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ABSTRACT 


This  report  contains  a  description  of  a  general  purpose,  6-degree 
of  freedom  terminal  homing  missile  simulation.  The  program  uses 
quaternions  for  the  coordinate  transformations  and  features  the  use 
of  subroutines  to  enter  seeker,  autopilot,  aerodynamic,  and  wind  models. 
The  program  is  written  in  FORTRAN. 
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Section  I  INTRODUCTION 


1.  Purpose  of  Proyam 

The  General  Purpose, 6-Degree  of  Freedom  Terminal  Homing 
Missile  reinitiation  Program  is  a  FORTRAN  program  designed  to  simulate 
the  dynamic**  of  terminal  homing  missile  systems  with  a  maximum  degree 
of  flexibility.  The  flexibility  is  achieved  through  the  use  of  user 
supplied  sunprograras  which  are  used  in  conjunction  with  the  main 
program. 

The  program  uses  quaternions,  as  opposed  to  Euler  angle  rotations, 
to  generate  the  coordinate  system  transformations.  The  quaternion 

approach  is  as  accurate  as  the  Euler  angle  approach^  and  has  the  advan¬ 
tage  of  avoiding  "gimbal  lock"  which  may  be  encountered  in  some  cases 
with  Euler  angle  rotations.  The  quaternions  do  not  need  FORTRAN  SIN 
or  COS  routines  in  the  computation  of  the  transformations  matrices  which 
results  in  a  potential  savings  in  the  computation  time. 


2.  Assumptions 

The  airframe  is  assumed  to  be  a  rigid  body  and  aeroelastic 
effects  are  not  included. 

The  airframe  is  assumed  to  have  a  plane  of  mass  syim'otry  coinciding 

with  the  vertical  plane  of  reference  (plane  defined  by  the  missile  x  and 

2  axes  in  Figure  1).  The  y-axis  is  therefore  a  principal  axis  and  the 

products  of  inertial  I  and  I  vanish.  Thus,  if  mass  asymmetries  are 

xy  yz 

to  be  simulated,  the  xy-plane  must  be  used  (x2  plane  must  be  a  plane 
of  symmetry). 

A  flat  nonrotating  earth  with  constant  gravity  is  assumed. 
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3.  Coordinate  Systems  Used  by  the  Frcgram 

The  program  uses  two  coordinate  systems.  The  earth  fixed 
(inertial  system)  and  the  missile  body-axis  system.  The  earth  fixed 
coordinate  system  is  assumed  to  be  fixed  to  a  flat  earth  with  the  z-axis 
of  a  right  hand  coordinate  system  pointing  down.  The  missile  body-axis 
system  is  a  right  hand  coordinate  system  with  the  x-axis  aligned  with 
the  longitudinal  axis  of  the  missile.  The  coordinate  system  is  fixed 
to  the  missile  and  rolls  with  it.  The  relationship  between  the  earth 


^Robinson,  A.  C.,  On  the  Use  of  Quaternions  in  Simulation  cf  Rigid 
Body  Motion,  WADC  TR  58-17,  December  1958. 
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Figure  1.  Relationship  between  the  Earth  and  Body-Axis  Systems 


fixed  and  body-axis  systems  is  given  in  Figure  1  in  terms  of  Euler 
angles.  However,  Euler  angles  are  not  used  in  the  computational 
structure  of  the  program  except  for  the  initial  conditions  and  when 
Euler  angle  printout  is  selected  as  part  or  a  standard  printout 
option  (Section  V). 

4.  Six  Degrees  of  Freedom  in  the  Simulation 

The  6  degrees  of  freedom  of  the  airframe  consist  of  three 
translations  along  che  inertial  X,  Y,  and  Z  axes  (earth  fixed  system) 
and  three  rotations  about  the  body- fixed  x,  y,  and  z  axes.  The  rota¬ 
tions  are  expressed  by  the  quaternion  e^  +  ie1  +  je£  +  ke^. 

(Appendix  A  contains  a  discussion  of  quaternions.) 


2 


Section  II  EQUATIONS  USED  BY  THE  PROGRAM 

1.  Variable  List 


The  following  is  a  list  of  variables  which  are  used  in  the 
equations  in  the  simulation  along  with  their  corresponding  FORTRAN 
program  variable  name. 


EQUATION  VARIABLE 

FORTRAN  VARIABLE 

VARIABLE  DEFINITION 

al 

A(l) 

Element  of  the  coordinate 
transformation  matrix 

a2 

A<2) 

Element  of  the  coordinate 
transformation  matrix 

a3 

A  (3) 

Element  of  the  coordinate 
transformation  matrix 

V# 

ALF4 

Angle  of  attack  in  the  pitch 
plane 

bl 

A  (4) 

Element  of  the  coordinate 
transformation  matrix 

b2 

A(5) 

Element  of  the  coordinate 
transformation  matrix 

b3 

A(6) 

Element  of  the  coordinate 
transformation  matrix 

P 

BETA 

Angle  of  attack  in  the  yaw 
plane 

C1 

A  (7) 

Element  of  the  coordinate 
transformation  matrix 

C2 

A(8) 

Element  of  the  coordinate 
transformation  matrix 

c3 

A(9) 

Element  of  the  coordinate 
transformation  matrix 

1 

DPI 

Wing  deflection  fin  1  (pitch) 

o 

Z 

DY2 

Wing  deflection  fin  2  (yaw) 

'  3 

DP3 

Wing  deflection  fin  3  (pitch) 

4 

DY4 

Wing  deflection  fin  4  (yawl 

D 

D 

Diameter  of  missile  body 

AX 

DELXE 

Missile  to  target  displace¬ 
ment  in  earth  X  direction 

DELYE 


Missile  to  target  emplace¬ 
ment  in  earth  Y  direction 


DEIZE 


DT 

DTMIN 


EMAX 

X(4) ,  EO 
X(5),  El 
X(6),  E2 
X(7),  E3 


EX 

FY 

FYE 

FZ 

FZE 

G 

HX 

HY 

HZ 

ITERA 

IX 

IY 

IZ 


Missile  to  target  displace¬ 
ment  in  earth  Z  direction 

Integration  step  size 

Minimum  allowed  integration 
step  size 

Maximum  allowed  error  in 
integration 

Ouaternion  parameter 

Ouaternion  parameter 

Quaternion  parameter 

Quaternion  parameter 

Constraint  error  [See  Equa¬ 
tion  (11-12)] 

Body-axis  component  of  force 
in  x  direction 

Body-axis  component  of  force 
in  y  direction 

Earth  system  component  of 
force  in  Y  direction 

Body-axis  component  of  force 
in  z  direction 

Earth  system  component  of 
force  in  Z  direction 

Acceleration  due  to  gravity 

Angular  momentum  of  missile 
about  body-axis  x-axis 

Angular  momentum  of  missile 
about  body-axis  y-axis 

Angular  momentum  of  missile 
about  body-axis  z-axis 

Number  of  integration  itera¬ 
tions 

Missile's  moment  of  inertia 
about  body-axis  x-axis 

Missile's  moment  of  inertia 
about  body-axis  y-axis 

Missile's  moment  of  inertia 
about  body-axis  z-axis 


I 


I 

2X 

IZX 

Missile's  product  of  inertia 
in  the  x-z  plane 

K 

K 

Arbitrary  constant  used  in 
quaternion  constraint 

- 

KE 

KE  =  (K)  (e) 

m 

MASS 

Mass  of  the  missile 

M 

MACH 

Mach  number  of  the  missile 

M 

X 

MX 

Moment  about  the  body-axis 
x-axis 

M 

y 

MY 

Moment  about  the  body-axis 
y-axis 

- 

MAXPT 

Maximum  allowed  number  of 
printouts 

M 

2 

MZ 

Moment  about  the  body-axis 
z-axis 

n 

s 

SX 

Number  of  seeker  and  auto¬ 
pilot  state  variables 

nt 

TX 

Number  of  target  state 
variables 

P 

X(l),  P 

Angular  rate  about  the  body- 
axis  x-axis 

- 

PRNTI 

Print  interval 

A 

PHI 

Euler  angle  rotation  phi 

- 

PS  I 

Euler  angle  rotatio.;  psi 

qs,  qsd 

QS,  GSR,  OSD 

Dynamic  pressure  terms 

q 

X(2),  Q 

Angular  rate  about  the  body- 
axis  y-axis 

_ 

RHO 

Air  density  constant 

r 

X(3) ,  R 

Angular  rate  about  the  body- 
axis  2 -ax  is 

R  . 
mm 

RMIN 

Miss  distance  (computer  after 
run  is  complete) 

- 

R(4) 

Integration  status  parameter 
(See  Appendix  B) 

- 

RM 

Range  to  the  target 

S 

S 

Reference  area  of  missile 

°y 

SIGY 

Line  of  sight  angle  in  earth 
coordinate  system  X-Z  plane 

5 


A 


A 


. .  . 


SIGZ 

SMAX 

X(I), 

14  5  I  s 

X(I),  ns 
and  I  £ 

TMAX 

THTA 

U 

X(8) 

UE 

V 

X(9) 

VE 

VM 

VS 

W 

X(10) 

WE 


Line  of  sight  angle  in  earth 
coordinate  system  X-Y  plane 

Maximum  error  in  integration 
computation 

Seeker  and  autopilot  state 
ns  +  13  variables 

+  14  s  I  target  state  variables 
is  +  nt  +  13 


Maximum  time  for  a  computer 
run  (real  time) 

Euler  rotation  Theta 

Body-axis  component  of  missile 
velocity  in  the  x  direction 

Earth  system  component  of 
missile  velocity  in  the  X 
direction 

Body  system  component  of 
missile  airspeed  in  the  x 
direction 

Body-axis  component  of 
missile  velocity  in  the  y 
direction 

i»rth  system  component  of 
r.<s>:ile  velocity  in  the  Y 

’»  ;  •.  ,.en 

E'.rtn  t  ’Stem  component  of 
mi.  silo  airspeed  in  the  Y 
ulrection 

Velocity  of  the  missile 

Velocity  of  sound 

Body-axis  component  of 
missile  in  the  z  direction 

Earth  system  component  of 
missile  velocity  in  the  Z 
direction 

Earth  system  component  of 
missile  airspeed  in  Z  direc¬ 
tion 


A 


•  C 


X(ll),  XE 


X (1 2^ ,  YE 


XrnK  ZE 


Earth  system  component  of 
wind  velocity  in  X  direction 

Earth  system  component  of 
wind  velocity  in  Y  direction 

Earth  system  component  of 
wind  velocity  in  Z  directior 

Body-axis  component  of  missile 
displacement  in  the  x  direc¬ 
tion 

Earth  system  component  of 
missile  displacement  in  the 
X  direction 

Earth  system  component  of 
missile  target  in  the  X 
direction 

Body-axis  component  of 
missile  displacement  ir.  y 
direction 

Earth  system  component  of 
missile  displacement  in  Y 
direction 

Karth  system  component  of 
target  displacement  in  Y 
direction 

Body -ax is  component  of 
missile  displacement  in  z 
direction 

Earth  system  component  of 
missile  displacement  in  Z 
direction 

Earth  system  component  of 
target  displacement  in  Z 
direction 


A  'dot"  over  a  variable  will  he  used  to  indicate  the  time  derivative 

of  the  variable.  For  example  X  is  the  time  rate  change  of  X.  The 
derivative  of  the  FORTRAN  variable  XH)  will  be  indicated  by  the  variable 
DX(T). 


2.  Initial  Conditions 

The  initial  conditions  for  the  simulation  are  given  in  terms 
of  the  equation  variables:  X,  Y.  Z,  u,  v,  w,  9,  :  ,  p,  q,  r,  Sj, 
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,  . . .  . 


t  f> 


S2,  •  •  .  Sn  «  Tj,  T2,  .  .  .  T^  wnich  are  supplied  by  the  user 
s  1 1 

equation  variables  to  be  integrated  are : 


Thus,  X,  Y,  Z,  eQ,  e^,  ancl  e^  must  be  computed  from  the  input  quanti¬ 
ties  supplied  by  the  user,  and  Equation  (A-33)  of  Appendix  A,  gives 

e_j  =  cos  l'/2  cos  0/2  cos  t>/2  -f  sin  -y/2  sin  9/2  sin  $/2 

e.,  =  cos  f/2  cos  9/2  sin  <t> / 2  -  sin  \i/2  sin  0/2  cos  $/2 

e2  =  cos  y/2  sin  0/2  cos  $/2  +  sin  ^/2  cos  0/2  sin  o/2 

e3  =  -cos  v/2  sin  9/2  sin  <t>/2  +  sin  v/2  cos  0/2  cos  o/2  .  (11-1) 


To  find  X,  Y,  Z  we  must  generate  a  transformation  matrix  which  will 

*  •  • 

take  the  velocities  u,  v,  and  w  (body-axis  system)  into  X,  Y,  Z  (inertial 
reference  system).  Thus,  using  Equation  (A-29)  of  Appendix  A, 


tUlU  Vli, 
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’&l’ 

hx(Sl"  \) 

&2 

MSi-  s2 . \) 

53 

h3(Si"  S2'-"  \) 

&4 

\(sv  s2,...snJ 

(iI-7) 

r  •  ~\ 

Ti 

fl(T1,  T?""’  \) 

T„ 

z 

= 

t2(tv  T2,...,  TnJ 

• 

T 

nt 

*JTl>  T2 . V) 

t  '  t/ 

(II-8) 

- 

XT 

3t(Tl '  T2 . \) 

= 

z2(t1’  \) 

ZT 

V"1  \) 

(II-9) 

e)  The  forces  and  moments  in  the  body  axis  system  -  f  ,  f 


S,  D,  M,  bjj  ^2>  ^3»  ^4’  ant*  P' 


fz>  My,  and  Mz  -  are  computed  by  user  supplied  equations  of  a,  £, 

The  moments  o:  inertia  I  and  I  ,  the 

x  y 

product  of  inertia  I  ,  and  the  mass  m  of  the  missile  are  also  supplied 
by  user. 

f)  Compute  the  following  equations  for  missile  angular 
accelerations: 


h  I  +  h 
X  z  z 


I  I  -  I 
X  z 


1 

zx 

2 

zx 
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h)  Compute  the  forces  in  the  inertial  system 


i)  Compute  the  accelerations  in  the  inertial  system 


U  =  F  /m 
x 

X  =  U 

V  =  F  /m 

V 

Y  =  Fz/m  +  g 

Z  =  W  .  (11-14) 

j)  Integrate  the  following  equations  variables: 


p,  q,  r,  U,  V,  W,  X,  Y,  Z,  eQ,  e1?  e2>  e3 

k)  Compute  the  components  of  wind  speed  W  ,  W  ,  W  in  the 

x  y  z 

inertial  system  by  user  supplied  equations. 

l)  Compute  the  inertial  components  of  missile  airspeed 


V  =  V  -  W 
e  y 

We  =  W  -  Wz  .  (11-15) 
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(Uj  V, 


where 


m)  Compute  the  body  axis  components;  of  missile  airspeed 
and  w). 


u 

aT  an  ao 

~  U  “ 

12  3 

e 

V 

b,  b,,  b. 

V 

12  3 

e 

w 

c,  c.  c_ 

W 

- 

1  2  3  _ 

e 

(11-16) 


2.2  2  2 

al  =  e0  +  ex  -  e2  -  e3 

*2  ‘  2(Cle2  +  V3) 

»3  -  2(ne3  -  v2) 

»1  '  2(V2  *  Vs) 


.  2,222 

b2  *  e0  +  c2  '  ‘l  -  e3 


b3  *  2(S2e3  +  e0ei) 

C1  *  2(ele3  +  V2) 

C2  *  2(V3  -  e0el) 

C3  =  e0  +  e3  *  el  *  e2  *  (II'17) 

m)  The  computational  loop  is  closed  by  going  back  to  Step  a. 
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Section  ill  FORTRAN  SUBROUTINES  SUPPLIED  BY  USER 


1.  General  Information 

The  user  must  supply  the  following  subroutines:  SEEKER, 

TARGET,  FOROM,  WRT,  STRPLT,  and  WIND.  Some  of  these  subroutines  may  not 
be  required  and  may  consist  of  only  a  DIMENSION,  a  RETURN,  and  an  END 
statement . 

Variables  not  in  the  argument  list  for  these  subroutines  must  be 
"transferred”  by  COMMON  statements.  The  main  program  contains  some 
standard  common  blocks  which  can  be  used  as  needed  (Paragraph  2  of 
Section  V) . 

2.  Subroutine  SEEKER 

Subroutine  SEEKER  is  used  to  model  the  dynamics  of  the  seexer 
and  autopilot  sections.  In  most  cases  common  block  ANG  will  be  required 
because  the  wing  deflections  will  be  needed  in  other  subroutines  and  in 
the  main  program  when  the  standard  printout  option  is  used.  Section  v 
of  this  report  contains  additional  information  on  the  standard  conmon 
blocks . 

Subroutine  SEEKER  should  contain  a  model  of  the  form  given  in  step 
d)  of  Paragraph  3  of  Section  II.  The  subroutine  should  define  the 
derivatives  of  the  state  variable  to  be  integrated  bv  using  the  FORTRAN 

variable  DX(I)  where  I  =  14,  15,  .  .  . ,  SX  +  13. 

For  example,  let  us  assume  that  the  seeker  and  autopilot  section 
are  defined  by  the  block  diagram  given  below  (the  s  is  a  Laplace  operator. 

Let  the  variables  shewn  in  Figure  2  be  defined  as  follows: 

Cy  -  Line  of  sight  angle  in  inertial  space  for  the  X-Z  plane  (rad) 

cz  -  Line  of  sight  angle  in  inertial  space  for  the  Y--Z  plane  (rad) 

*p  -  Pitch  wing  command  (rad) 

r-  -  Yaw  wing  command  (rad) 

K  -  Seeker  gain  (rad /sec /rad) 

Kn  -  Navigation  gain. 
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The  corresponding  FORTRAN  program  is: 


SUBROUTINE  SEEKEr<TImE*X,DX) 

REAL  INPUT 
REAL  KE 
REAL  K  »KN 

DIMENSION  X(i )»0X(1 ) 

COMMON  /ANG/ALFA,BCTA.SIGY,SIGZ,DP:WDY2.DP3»DY4 
KN--8.5 

DX(14)=K*(SIGY-X(14>) 

IX(1&)=K#(SIGZ-X(15)  ) 

DP1=KN*DX (14) 

DY2=KN*DX(1?5) 

IF ( AgS(DPl) *GT, .1745)  DP1=SIGN( .174b,0Pl) 
IF(ABS(DY2)  .GT.  .1745)  1)Y2=SIGN<  .1745.DY2) 
IFITIME.LT.50,)  OH1=DY2=0. 

DY4=pY2 
DP  5-  DP  1 

RETURN 

END 


Note  that  the  wing  commands  DPI,  DP3,  DY2,  and  DY4  may  be  defined  in 
any  manner  the  user  wishes,  however,  the  user  must  be  consistent  with 
the  wing  command  defined  in  the  user  supplied  FOROM  subroutine.  Note 
also  that  the  FORTRAN  variable  SX  is  equal  to  two  because  there  are  two 
seeker  states  [X(14>  and  X(15)]. 


3.  Subroutine  TARGET 

Subroutine  TARGET  is  used  to  model  the  target  dynamics.  The 
model  should  have  the  form  of  the  model  given  in  step  d).  Paragraph  3 
of  Section  II.  The  subroutine  es  a  general  rule  should  contain  the 
standard  common  block  DISPL  (Section  V).  The  derivatives  of  the  state 
variables  to  be  integrated  must  be  defined  by  the  FORTRAN  DX(I),  where 
I  =  13  +  SX  +  1,  13  +  SX  +  2,  ...,  13  +  SX  +  TX.  The  target  locations 
in  earth  coordinate  system  (XTE,  YTE,  and  ZTE)  must  also  be  defined. 

Subroutine  TARGET  can  be  illustrated  by  a  simple  example  where  the 
target  is  located  at  X  =  46,000  ft,  Y  =  500  ft,  Z  =  0  ft  in  the  earth 
coordinate  system  and  is  moving  with  a  velocity  of  5  ft/sec  in  the  Y 
direction. 


SUBROUTINE  TaRGcT( TlME,X,DX) 
dimension  x<i>,nx<i> 

COMMON  /BISPl/  U,V,W,DELXE,nELYE.DELZt»XTE,vTt,ZrE,BM 
DX(16)=5. 

XTE  s4  6000 , 

YTEs500,*X(16> 

ZTEs0. 

return 

END 

Note  that  there  's  one  state  variable  in  the  subroutine  and  the  FORTRAN 
variable  TX  is  ecual  to  1. 


4.  Subroutine  FOROM 

The  purpose  of  subroutine  FOROM  is  to  supply  the  forces, 
moments,  moments  of  inertia,  product  of  inertia,  and  mass  to  the  main 
program.  The  forces  are  in  the  body-axis  system. 

The  positive  conventions  for  forces,  moments,  and  angular  rates  for 
the  body-axis  system  are  shown  in  Figure  3. 


Figure  3.  Positive  Conventions  in  the  Body-Axis  System 

Subroutine  FOROM  must  contain  common  blocks:  FRMO,  AIR,  ANG,  and 
0TKER2  for  most  applications  (Section  V)  likely  to  be  encountered. 

The  use  of  subroutine  FOROM.  can  be  illustrated  by  a  simple  example: 


SUBROUTINE  FOROM(TIME.x.DX) 

INTEGER  SX.TX 
REAL  KE.K 

REAL  MACH»MX,mY,MZ,IX,  IY,IZ.  IZX.MASS 
DIMENSION  X<l)tDX(l) 

COMMON  /  ANG/  ALpA  ,BET  A#  SI  GY  rS  IGZ»  DPI » RY2#  DP3»  DY4 
COMMON  /0THER2/  THTA ,RSl ,PHl ,KE* S#D, SX, TX»K , G» MAXPT 
COMMON  /AIR/  MACm, vm,vs»qs»osd.qsr,rho 
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COMMON/FRMO/  ?X » FY , FZ . MX , MY , MZ , I X , I Y , I Z , I ZX . MASS 

IY=IZ=4.77 

IX-.202 

MASS=4 .671 

IZX-0. 

Ss  .1 96  5 
P=.5 
CD  Os  .2  5 
CMA=-2  3. 

CN  A=  15  . 

CMD=  18  . 

CMOs-263. 

CALL  QSDSUB 
FX=- ,5*QS*CDO 
FYs-,5*QS*CNA*tiETA 
FZ=-,5«QS*CNa«ALFA 
QsX(2) 

RSX (3 ) 

HX=?. 

mYs.5*QSd*(CMA*AlFa+CMD*DP1*D«CMQ*Q/(2 ,*VM > ) 

MZ  X.  5«  QS  D*  (“  C”  A*  PE  TA  +C  MD  *D  Y2  +D  #C  MO  «P  /(  2.  *V  M)  ) 

RE  TURN 
END 

Note  in  the  example  the  aerodynamic  coefficients  are  assumed  to  be 
constants.  In  most  simulations  this  assumption  would  not  be  valid  and 
aero  data  must  be  stored  in  tables  as  a  function  of  Mach  number.  For 
some  tail  controlled  missiles  where  the  angles  of  attack  are  likely  to 
be  relatively  large,  it  may  be  necessary  to  compute  the  aerodynamic 
coefficients  as  a  function  of  three  independent  variables  -  Mach  number, 
angle  of  attack,  and  wing  deflection.  The  problem  is  basically  one  of 
interpolating  to  find;  an  accurate  aerodynamic  coefficient  for  any  given 
set  of  independent  variables.  To  help  solve  this  problem,  an  interpola¬ 
tion  routine  has  been  provided  as  part  of  the  main  progran  which  will 
handle  one,  two,  or  three  independent  variables.  A  description  of  the 
interpolation  routine  can  be  found  in  Section  IV.  Subroutine  QSDSUB 
is  a  subroutine  supplied  by  the  main  program  which  computes  the  dynamic 
pressure  terms  (Section  IV). 

5,  Subroutine  WRT 

The  purpose  of  subroutine  WRT  is  to  provide  a  way  to  output 
certain  variables  not  in  the  standard  printout.  As  will  be  discussed 
in  Section  V,  there  is  an  input  data  card  which  provides  the  user  with 
the  following  output  options:  standard  printout  only,  printout  provided 
by  the  user  in  WRT  only,  or  the  printout  provided  by  the  standard  print¬ 
out  plus  the  output  generated  by  WRT.  Subroutine  WRT  must  as  a  minimum 
contain 
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SUBROUTINE  Wm(TlM£,X,DX) 

dimension  x(i).Dxd) 

RETURN 

END 

For  this  exanple,  the  standard  printout  option  would  normally  be  used. 
All  variables  to  be  printed  out  by  WRT  with  the  exception  of  the  arrays 
X  and  DX,  and  the  variable  TIME  must  be  transferred  to  WRT  by  COMMON 
statements. 


6.  Subroutine  STRPLT 

The  purpose  of  the  user  supplied  subroutine  STRPLT  is  to 
store  output  data  in  arrays  as  required  and  to  provide  enough  flexibility 
so  that  plots  can  be  made  with  either  the  line  printer  or  with  some 
other  output  device  such  as  the  Tektronix  4002A  graphics  terminal. 

A  line  printer  plot  will  be  generated  by  the  following  example 
program.  Subroutine  PLOT  is  discussed  in  Section  IV.  The  line  printer 
plot  generated  by  this  program  is  given  in  Section  VI. 


subroutine  s tr pl t (  time.x.pa,  mss.  IpL3T; 

DIMENSION  X(  1)  »DX(  l) 

DIMENSION  A(iwflZ) ,0(500) 

CO  “K  ON  /0  TH  ER  1/  I  Th  RA  ,D  T,  DT  r-j,  FM  AX  ,  \  »S  MA  X,  TM  AX  ,P  RN  TI 
FORMAT (lKl#5dX,*ZE  V S  TlMF«//  ) 

IF  {  ITERA  .Nb.  1)  GO  TO  1 
1  =  4* 

TS*3. 

STFP«1. 

CO  NT  IiJUE 

I F( TJME . L  r . T$ , AND . IPLOT  «EQ . 0)  00  TO  3 

ts=t  i«e*sTgp 
l  =  i*i 
A( I)*TIME 
3< I )sX(li) 

IF  (I  PLOT  ,Fq,  j)  30  TO  3 

WR I TE ( 6»  4 ) 

Ns  I 

DO  5  Jsl.N 
A  (  J-M>s3{  J) 

CALL  PLOT (A  »N  »2 .2  *  2  ) 

RETURN 

END 
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The  FORTRAN  variable  IPLOT  is  set  to  zero  during  the  flight  phase  of  the 
simulation  and  is  set  to  one  when  the  run  is  complete  and  plotting  can 
commence.  IMISS  =  0  when  the  miss  distance  is  less  than  50  feet. 

IMISS  =  1  when  the  miss  distance  is  greater  than  50  feet  (RMIN  has  no 
meaning) . 


7.  Subroutine  WIND 

The  purpose  of  the  user  supplied  subroutine  WIND  is  to  provide 
a  method  for  placing  a  wind  model  in  the  simulation  when  it  is  required. 

Subroutine  WIND  should  contain  a  minimum  of  the  following 

SUB^OUTI  NE  W  IND<  Tl  NE  ,X  ,DX»  U#  V»  V  > 

DIMENSION  X<1>  ,DX(1> 

RETURN 

END 

The  variable  U,  V,  and  W  are  components  of  missile  velocity  in  the 
earth  coordinate  system.  So  to  enter  a  wind  model;  U,  V,  and  W  must 
be  "replaced"  with  inertial  components  of  missile  airspeed.  For  example, 
if  the  wind  is  blowing  in  the  inertial  Y  direction  with  a  velocity  of 
10  ft/sec,  the  subroutine  would  take  the  following  form: 

subroutine  wind  <ti**e.x,dx,u,v,w> 
dimension  x<i>  ,qx(i) 

V=V-  IP*. 

RcT-JRN 

end 
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Section  IV  FORTRAN  SUBROUTINES  SUPPLIED  BY  PROGRAM 


1.  General  Information 

There  are  certain  subroutines  available  with  the  main  program 
which  may  be  called  by  the  user  supplied  subroutines.  The  subroutines 
simplify  the  task  of  programming  in  that  a  group  of  statements  may  be 
replaced  by  a  single  call  statement  in  the  user's  program.  The  following 
sections  define  the  function  of  each  subroutine  which  can  be  called  by 
the  user. 


2.  Subroutine  LIM 

Subroutine  LIM  contains  a  limiter  model  and  has  a  calling 
statement  of  the  form 


CALL  LIM  (INPUT,  OUTPUT,  A,  B) 


where 


OUTPUT  =  INPUT  when  A  £  INPUT  £  B 
OUTPUT  =  B  when  INPUT  >  B 
OUTPUT  =  A  when  INPUT  <  A 


(IV-l) 


NOTE:  INPUT  is  a  real  variable  in  subroutine  LIM. 

3.  Subroutine  LIMSTA 

Subroutine  LIMSTA  contains  a  limiter  model  which  modifies  the 
input  quantity  as  well  as  limits  the  output.  The  calling  statement  is 
of  the  form 

CALL  LIMSTA  (INTUT,  OUTPUT,  A,  B) 

where 


OUTPUT  =  INPUT  when  A  £  INPUT  £  B 
OUTPUT  =  INPUT  =  B  when  INPUT  >  B 

OUTPUT  =  INPUT  =  A  when  INPUT  <  A 

INPUT  is  a  real  variable  in  subroutine  LIMSTA. 


(IV-2) 


4.  Subroutine  OETEC 


Subroutine  DETEC  contains  a  simplified  detector  model  with  a 
characteristic  as  shown  in  Figure  4. 


Figure  4.  Model  of  Detector  Characteristic 
The  calling  statement  is  of  the  form 

CALL  DETEC  (INPUT,  OUTPUT,  A,  B,  C,  D) 

where  the  arguments  are  defined  in  Figure  4. 

NOTE:  INPUT  is  a  real  variable  in  subroutine  DETEC. 


5.  Subroutine  DEADSP 

Subroutine  DEADSP  contains  a  medal  of  an  element  with  dead 
space  as  shown  in  Figure  5.  The  calling  statement  is  of  the  form 

CALL  DEADSP  (INPUT,  OUTPUT,  A,  B) 


where  the  arguments  are  defined  in  Figure  5. 

NOTE:  INPUT  is  a  real  variable  in  subroutine  DEADSP. 


S.  Subroutine  LAG 


transfer 


Subroutine  IAG  contains  a  model  of  a  first  order 
function  of  the  form 


leg  with  a 
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(IV -3) 


OUTPUT (S)  _  1 

INPUT(S)  ~  ts  +  1 


Figure  5.  Dead  Space  Model 
The  calling  statement  is  of  the  form 

CALL  LAG  (INPUT,  OUTPUT,  X,  DX,  INDF  T) 

where  X  is  the  array  of  state  variables  used  by  the  main  program 
(similarly  DX  is  the  array  consisting  of  the  derivative  of  X),  INDEX 
xs  an  integer  defining  tne  element  of  X  which  defines  the  state  of  the 
first  order  lag,  and  T  is  the  time  constant  t  defined  in  Equation  (IV-3) . 
Note  that  a  first  order  lag  model  requires  only  one  state  variable  and 
that  INPUT  is  a  real  variable  in  subroutine  LAG. 


7.  Subroutine  SECORD 

Subroutine  SECORD  contr.ins  a  model  of  a  second  order  linear 
system  with  a  transfer  function  of  the  form 


OUTPUT (S) 
INPLT(S) 


w_ 


n 


S2  +  2*u  S  +  „2 
n  n 


(IV-4) 


The  calling  statement  is  of  the  form  | 

i 

t 

CALL  SECORD  (INPUT.  OUTPUT,  X,  DX,  INDEX,  ETA,  WN)  \ 


A 

1 

1 


$ 
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where  ETA  =  |  and  WN  =  un  [Equation  (IV-4)]  and  INDEX  is  an  integer 

defining  the  element  of  the  array  X  (Section  IV)  which  corresponds  to 
the  first  state  variable  of  the  model  for  the  second  order  system.  Note 
that  two  states  are  required  to  model  a  second  order  system  and  that 
INPUT  is  a  real  variable  in  subroutine  SECORD. 


&  Subroutine  LDLAG 

Subroutine  LDLAG  contains  a  model  of  a  standard  lead-lag 
compensation  network  with  a  transfer  function  of  the  fora: 

OUTPUT (S)  T1S  *  1 

INPUT (S)  =  t2s  +  1  .  (IV-5) 


The  calling  statement  is  of  the  form 

CALL  LDLAG  (INPUT,  OUTPUT,  X,  DX,  INDEX,  Tl,  T2) 

where  INPUT,  OUTPUT,  X,  DX,  and  INDEX  art  defined  similarly  to  the 
definitions  given  ir:  Section  IV  and  where  Tl  =  and  T2  «  t2  [Equation 

(IV-5)].  Note  that  one  state  variable  is  required  to  model  a  lead-lag 
network  and  that  INPUT  is  a  real  variable  in  subroutine  LDLAG. 

9.  Subroutine  GYR02 

Subroutine  GYR02  contains  a  model  of  a  seeker  gyro  which  can 
be  torqued.  The  calling  statement  is  of  the  form 

CALL  GYR02  (TGY,  TGZ ,  X,  DX,  INDEX,  IX,  ITP,  ITT,  WS)  . 


Consider  Figure  6. 


TGY 


x.y,i  •  BODY-AXIS,  AXES 

X_Y„  Z.  -  GYRO  COORDINATE 

w  9.  • 

NOTE:  Y#  AND  Z#  AXES  ARE 


|  TGZ 

Coordinate  System 


SYSTEM  AXES 
NOT  SHOWN 
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Let  x,  y,  and  z  be  the  coordinates  of  the  body-axis  system  (fixed  on  the 
missile)  and  let  6  and  \}r  define  the  location  of  the  gyro  spin  axis.  Then 
the  arguments  of  subroutine  GYR02  are  defined  as  follows: 

TGY  -  Torque  on  the  gyro  as  defined  in  Figure  6 
TGZ  -  Torque  on  the  gyro  as  defined  in  Figure  6 

X  -  Array  of  state  variables  [X  (INDEX)  =  9  and  X (INDEX  +  1)  =  \|r] 
IX  -  Axial  moment  of  inertia  of  the  gyro 
ITP  -  Transverse  moment  of  inertia  in  the  pitch  plane 
ITy  -  Transverse  moment  of  inertia  in  the  yaw  plane 
WS  -  Spin  rate  of  the  gyro 

INDEX  -  An  integer  defined  as  in  Section  IV. 

This  gyro  model  does  not  model  nutation  of  the  gyro  and  thus  is 
suitable  for  digital  integration  using  a  relatively  large  step  size. 

Note  that  two  states  are  required  to  model  the  gyro.  The  torques  TGT 
and  TGZ  are  transformed  to  torques  in  the  gyro  coordinate  system  which 
has  an  X-axis  aligned  with  the  gyro  spin  axis  and  which  does  not  rotate 

O 

or  spin  with  the  gyro  wheel.  The  equations  used  in  the  model  can  be 
obtained  from  the  listing  of  subroutine  GYR02  given  in  Appendix  C. 

10.  Subroutine  PLOT 

Subroutine  PLOT  was  intended  to  be  called  in  subroutine 
STRPLT  (Section  III)  and  can  be  used  to  plot  out  information  using  the 
line  printer.  Section  VI  contains  an  example  of  a  trajectory  (altitude 
versus  time)  made  with  subroutine  PLOT. 

The  calling  statement  is  of  the  form 

CALL  PLOT  (A,  N,  M,  NL,  NS) 


where 


A  -  Matrix  of  data  to  be  plotted.  The  "irst  column  represents 
the  base  variable  (e.g.,  time),  and  the  successive  columns 
are  the  cross  variables  (the  maximum  number  of  cross  variables 
is  nine). 

N  -  Number  of  rows  in  A 
M  -  Number  of  columns  in  A 

NL  -  Number  of  lines  in  plot.  If  NL  =  0  specified,  50  lines  are 
used  (one  standard  line  printer  page) 
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NS  -  Code  for  sorting  the  base  variable  into  ascending  order. 

If  NS  =  0,  the  base  variable  is  in  ascending  order  and  sorting 
is  not  required.  If  NS  =  1  sorting  is  necessary  and  the  base 
variable  will  be  placed  in  ascending  order. 


1 1 .  Subroutine  TABLE 

Subroutine  TABLE  is  an  interpolation  subroutine  which  will 
handle  functions  with  one,  two,  or  three  independent  variables. 

The  call  statement  is  of  the  form 

CALLTABIE  (N,  ANSWER,  WT,  X,  XT,  NX,  NPX,  Y,  YT,  NY,  NPY ,  Z,  ZT,  IE,  NPZ) 


where 


N  - 
ANSWER  - 
WT  - 
X  - 
XT  - 
NX  - 
NPX  - 
Y  - 
YT  - 
NY  - 
NYX  - 
Z  - 
NZ  - 
NPZ  - 


The  number  of  independent  variables 
Dependent  variable  corresponding  to  (X,  Y,  Z) 

Table  of  dependent  variables  corresponding  to  (XT,  YT,  ZT) 

Independent  variable  X 

Table  of  independent  X  values 

Number  of  points  in  XT 

Number  of  points  to  be  used  for  X  interpolation 

Independent  variable  Y 

Table  of  independent  Y  values 

Number  of  points  in  YT 

Number  of  points  to  be  used  for  Y  interpolation 
Independent  variable  Z 
Number  of  points  in  ZT 

Number  t  f  points  to  be  used  for  Z  interpolation. 


The  use  of  subroutine  TABLE  can  be  illustrated  by  a  simple  example. 
Assume  that  Table  I  is  a  table  of  values  for  the  axial  drag  coefficient 


The  data  can  be  stored  in  tables  by  means  of  a  data  statements 


DATA  ((  (CD0T<  I*  J»  K)  » I  =1  «b  ),  J=l,  4,  K=l»  2)  / 
X , 8 , . 7, .6, .4, .3, 

X  ,6  » .  5»  .4, .3*  .2  » 

X,5, . 4 » ,  2  #  . 2 » . 2 > 
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X  ,  4  , ,  4  > .3, .2,  .2, 

X.  7,  ,6  # .  5»  .4,. 3, 

X ,5  » ,  4, .3 , ,  2*  .2# 

X3.»3,#3.»2,*2** 

X, 3, ,3* .3, .4, .5/ 

DATA  DEL< I >  *  I =1 »4 )/-l5, »-10 « » -5 . , 0. / 

DATA  (MACH (I  ) » Isl»2)/,4» .8/ 

DATA  ULFATCI  >,  Isl,5>/0.,5,  ,10.  ,15.  ,20,/ 


Table  I.  Data  Table  for  Drag  Coefficient  C^ 


0 

5 

Angle  of  Attack 
(deg) 

10 

15 

20 

Wing 

Deflection 

(deg) 

Mach 

No. 

- - 

0.8 

0.7 

0.6 

0.4 

0.3 

-15 

0.4 

0.6 

0.5 

0.4 

0.3 

0.2 

-10 

0.4 

0.5 

0.4 

C.2 

0.2 

0.2 

-5 

0.4 

0.4 

0.4 

0.3 

0.2 

0.2 

0 

0.4 

0.7 

0.6 

0.5 

0.4 

0.3 

-15 

0.8 

0.5 

0.4 

0.3 

0.2 

0.2 

-10 

0.8 

0.3 

0.3 

0.3 

0.2 

0.2 

-5 

0.8 

0.3 

0.3 

0.3 

0.4 

0.5 

0 

0.8 

The  program  would  contain  a  call  statement  of  the  form 


CALL  TA8Lt<3,CD0,CD0T,  ALFA  ,  A(.F  AT  ,  5,2,  DPI  ,DEu ,  4 ,2  ,  MACH,  *T  ,  2 , 2  ) 

Thus  given  ALFA,  DPI,  and  MACH,  CDO  will  be  computed  by  interrelation. 
Note  that  NPX,  NPY,  and  NPZ  are  equal  to  two.  This  implies  chat  all 
interpolation  will  be  linear.  Note  also  that  the  data  tables  ere 
arranged  in  ascending  order  of  magnitude  for  the  independent  »<.riaojes 
and  that  the  angle  of  attach  table  values  correspond  to  the  !  subscript 
of  the  data  statement  for  CDOT(I,  J,  K) .  The  J  and  K  subscripts  corre¬ 
spond  to  the  wing  deflection  and  Mach  number  table  values  respectively. 
When  N  =  2,  the  variables  Z,  ZT,  IC ,  N VZ  are  all  auasy  variables.  When 
N  =  1,  the  variables  Y,  YT,  NY,  NYZ ,  Z,  ZT,  NZ,  NPZ  are  all  dummy 
variables. 
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12.  Subroutine  QSDSUB 


Subroutine  QSDSUB  contains  computations  for  the  standard 
aerodynamic  pressure  terms.  The  calling  statement  is  of  the  form 

CALL  QSDSUB 

The  cotonon  block  AIR  must  be  included  in  the  subroutine  which  call 
QSDSUB.  The  variables  in  common  block  AIR  are  computed  as  follows 


QS  -  o  S 
QSD  =  p  V2  S  D 

Xu 

QSR  =  p  V  S  D2 
w  m 


Section  V  INPUT  AND  OUTPUT  DATA 


1.  Input  Data 


The  following  is  a  list  of  input  data  cards  required  by  the 

program: 


a)  SX  TX  (read  in  with  215  format) 

b)  TAG(I) ,  I  =  1,  8  (read  in  with  8A10  format) 

c)  IOPTN,  INTOPT  (read  in  with  215  format) 

d)  TIME,  PHI,  THTA,  PSI  (read  in  with  8F10.0  format) 

e)  U,  V,  W,  X(l),  X(2) ,  X(3)  (read  in  with  8F10.0  format) 

f)  X(ll),  X(12) ,  X(13)  (read  in  with  8F10.0  format) 

g)  X(I) ,  I  =  14,  13  +  SX  (read  in  with  8F10.0  format) 

h)  X(I),  I  =  13  -f  SX  +  1,  SX  +  TX  +  13  (read  in  with  8F10.0 

format) 

i)  DT,  OTMIN,  F.M£X,  PRNTI,  TMAX  (read  in  with  8F10.0  format) 

j)  MAXPT  (read  in  with  215  format) 

k)  (Blank  card) 


where 


SX  is  the  number  of  state  variables  in  subroutine  SEEKER 

TX  is  the  number  of  state  variables  in  subroutine  TARGET 

TAG(I),  1  =  1,  8  is  a  80  column  title  for  output  (may  be  a  blank  card) 

IOPTN  is  the  print  option  (if  IOPTN  =  0  standard  printout;  if 
IOPTN  =  1,  printout  by  subroutine  WRT;  and  if  IOPTN  =  2, 
standard  printout  plus  printout  provided  by  WRT) 

INTOPT  is  the  integration  routine  option  (if  INTOPT  =  0,  Runga- 
Kutta  fourth  order  used  for  integration;  if  INTOPT  =  1, 
variable  step  Runga-Kutta  Merson  is  used  for  integration; 
if  INTOPT  =  2,  Hamming  Predictor  Corrector  is  used  for 
integration) .  See  Appendix  B  for  discussion  of  integration 
routines. 

TIME  is  the  initial  value  of  time 

PHI  is  the  initial  value  of  the  Euler  angle  <t>  (Figure  1) 

THTA  is  the  initial  value  of  the  Euler  angle  9  (Figure  1) 

PSI  is  the  initial  value  of  the  Euler  angle  \{r  (Figure  1) 

U  is  the  initial  value  of  the  missile  velocity  in  the  body-axis 
x-direction 
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V  is  the  initial  value  of  the  missile  velocity  in  the  body-axis 
y-direction 

W  is  the  initial  value  of  the  missile  velocity  in  the  body-axis 
z-direction 

X(l)  is  the  initial  value  of  the  body  rate  p  about  the  bi.dy-axis 
x-axis 

X(2)  is  the  initial  value  of  the  body  rate  q  about  the  body-axis 
y-axis 

X(3)  is  the  initial  value  of  the  body  rate  r  about  the  body-axis 
z-axis 

X(ll)  is  the  initial  displacement  of  the  missile  in  the  earth 
coordinate  system  X-direction 

X(12)  is  the  initial  displacement  of  the  missile  in  the  earth 
coordinate  system  Y-direction 

X(13)  is  the  initial  displacement  of  the  missile  in  the  earth 
coordinate  system  Z-direction 

X(I),  I  =  14,  13  +  SX  are  the  initial  values  for  the  state  variables 

in  subroutine  SEEKER 

X(T),  I  =  13  +  SX  +  1,  SX  +  TX  +  13  are  the  initial  values  for  the 

state  variables  in  subroutine 
TARGET 

DT  is  the  step  size  to  be  used  for  integration 

DTMIN  is  the  minimum  step  size  to  be  used  for  the  variable  step 
size  integration  methods.  The  constant  has  no  significance 
for  the  fixed  step  Runga-Kutta  integration  routine  but  DTMIN 
must  still  be  defined. 

EMAX  is  the  maximum  error  to  occur  in  the  integration  routines. 

The  step  size  will  be  reduced  until  the  error  is  less  than 
EMAX  or  the  minimum  step  has  been  reached  (DTMIN)  (DTMIN  used 
for  Runga-Kutta-Merson  only) 

PRNTI  is  the  print  interval  in  seconds  (PRNTI  must  be  greater  than 
or  equal  to  DT) 

TMAX  is  the  minimum  value  for  the  FORTRAN  variable  TIME.  The  run 
is  terminated  and  the  next  data  set  is  read  in  if  TIME  >  TMAX. 

MAXPT  is  the  maximum  allowable  number  of  printouts.  After  the  run 
is  complete  for  one  data  set,  the  program  is  initialized, 
and  a  new  data  set  is  read  in.  Thus,  data  sets  may  be 
"stacked"  in  sets.  The  program  can  be  terminated  by  setting 
SX  =  TX  =  0.  This  can  be  accomplished  by  placing  a  blank 
card  after  the  last  data  set. 
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2.  Output  Data 


a.  Standard  Common  Blocks 


The  standard  common  blocks  which  are  part  of  the  main 

program  are 

CCFMCN  /TRANS/  A(9) 

CCKMCN  /OISPL/  U,  VfW DELXE,  DEI  YE  .TEL  ZE  ♦  XTE.YTE  ,ZTE  ,RH 
CCRMCN  /ANG/ALFA,  oft A,SIGY,SIG7» CPi, DY2  ,0P3» OY4 
CCKMCN  /AIR/  HACH,VH,VS,QS,OSD,QSR,RHC 
CCtMCN  /FPMO  /  FX»FY  ♦  FZ»MX*MY»MZ»IX?  IZ»IV  »IZX»  HASS 
CCFMCN  /0THEF1/IT ERA ♦DT»DTMINtEMAX,SMAX»T^AX  *PRNTI 
CCMM0N/0THER2/  TH  TA,  PSI  ,PHI  ,  KE  ,  S  ,D  ,S  X  ,  TX  ♦  K  ,G  *  HAXPT 
CCKMCN  /1NT/  RCA) 

COMMON  /MIS/  RMIN 


These  common  blocks  may  be  used  to  transfer  variables  from  the 
main  program  to  various  user  supplied  subprograms  (and  vice  versa)  and 
from  user  supplied  subroutine  to  user  supplied  subroutine.  The  variables 
given  in  the  COMMON  statements  are  defined  in  Section  II. 

The  following  variables  must  be  declared  to  be  real  when  a  common 
block  containing  them  is  used:  MX,  MY,  MZ ,  IX,  IY,  IZX,  MASS,  KE,  and 
K.  The  variables  SX  and  TX  must  be  declared  to  be  integer  variables 
when  common  block  0THER2  is  used. 

b.  Output  Options 

The  opt  Lons  available  for  outputing  have  been  cover'd  in 
Section  II.  A  sample  of  the  standard  printout  is  given  in  Section  VI. 
Note  that  Euler  angles  are  provided  in  the  standard  printout.  These 
angles  are  computed  only  at  the  print  interval  from  the  quaternions 
(Appendix  A) . 


Section  V!  EXAMPLE  PROGRAM  AND  OUTPUT 


1.  Subprograms  Supplied  by  User 

The  following  listing  is  comprised  of  the  listing  front  examples 
used  in  the  previous  sections . 


SUBROUTINE  SEEKER(TIME,XfO**  . 

REAL  INPUT 
REAL  Kt 
REAL  K,<N 

DIMENSION  X(l>  »  GX  ( 1) 

COMMON  /ANG/4LFA,BETA,SIGY,SIG7,  0r'l,aY2»DP3«DYU 
KN=3 . 5 
K  =  ld  . 

0X(l‘t)=<MSIGY-X(14)  ) 
r«X(i5)=<*(SIGZ-X(  15)  ) 
nPl=<N*GV (14) 
nV2=KN»GX (15) 

Ic  (A BS (QP1 )  •  GT  ..1745)  0P1=3IGN(  .17V5,0PJ  ) 

IF(ARS(3VE>.G1'  ..1745)  13  Y2=S  I GNC  .  1745  ,  DY2 ) 

TF(TI^fc.LT.53.  )  nPl=OY2  =  0. 

nVi.=  GY2 

GP3=QP1 

RETURN 

ENG 


SUBROUTINE  cOR  ON ( TINE  »  X , GX) 

INTE GE°  SX.TX 
REAL  <t,K 

RFAL  MA3H,KX*MY,MZ,IXtIY,IZ»IZX, MASS 
DIMENSION  Ml)  (l) 

COMMON  /AN&/ALFA,8ETA,StGY,SIGZ,  0P1 tOY2 1 t>P3,0Y4 
COMMON  /OTHER?/  TH7A , PSI , PMI , KE, S ,DtSX ,T X *K»G tMAXPT 
COMMON  / AIE/  MACH,VM,VS,QS,aSn,3SR,RHU 
COmMON/cRMO/  FX,FY.FZiMX,MY,M7»IX ,  IY  ,  IZ  ,  I ZX  ,  M  ASS 
IY  =  IZ=4.  77 
TX=. 202 
MASS -4.6  71 
IZX=d  . 

S= . 1 95B 
0=  .5 
CDO= .25 
CMA=  -23. 

CNA=15. 
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C*0=18. 

CM0=-283. 

call  osasun 

FX=-.5*3S*C30 
FY=- , 5*3S*CNA*PETA 
FZ=-  .  O^QS+ONA*  ALFA 
Q=X( 2) 

R=X( 3) 

MX-=3  . 

MY=.5*QSnMCMA*ALFA+CMO*OPl  +  D*CMQ*G/(2.*VM>  J 
H7=, 5+QSO* (-CM A*RETA+CMD*nV2+D*CMQ*R/ <2.*VM> ) 
RfT'JpM 
END 


SiiBRO'JT I  NF  STRPLT  (TIMFiXtB**  IMIS  S  * IPL0TI - 
Ot^NSION  X(  1)  ,0X  (1) 

3IHENSI3N  A(lJu3>  ,R(5C;> 

C0MM0N/3Th£R1/  ITf.RA,DY,DTMlN,f  MAX,N,SMAX,TMAX,PRNTI 
A  FORMAT (IH1 » 5 JX »*ZE  VS  TIMF*//  > 

I F  ( I  T cPA  «  NE  .  1)  GO  TO  1 
I  =  u 
T  S—  0  . 

STE3  =  1.. 

1  -CONTI  4UE 

IF<TH£.LT.T3«  AND  •  I°LOT.  trQ.  •) >  GO  To  3 

ts=time+step 
1=1+1 
A  { I)  =  TIME 
R(I)  =  <<1  J> 

IF(IPLCT.EP.  0)  GO  TO  3 
WRITE  {‘  ,hJ 
N=  I 

00  5  J  =  1  »  N 
6  A(J+N)=3(J> 

CALL  PLOT ( A ,  N,  2, 3  ,  G) 

3  RETURN 

FNC 


.e 


ll 

s 


fH 


SUSRO0TI  ME  TARGET  (TIME, X,DX> 

OKIE MSION  X{  1)  ,QX  (!) 

COMMON  /QISPL/  U , V,W , OELXE , Oe  LYE  , OELZE ,XT£,YTE»ZTE,RM 
OX (1 6) =5  » 

XTE=46GOO  . 

YTE=5  3u.  +  X  (16) 

Z  T  F  =  •» , 

RETURN 

END 
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SUBROUTINE  wrutime,  x,ox> 

DIMENSION  Xfl)  ,DX  C.U 

pETJ.Rn 

END 


SUBROUTINE  <HP(TIME  ,X  ,OX,'J,V  ,W> 
0IHENSI3N  XC  i)  ,0X  (1) 

RETURN 

ENG 

2.  Input  Data 


The  following  initial  conditions  will  be  assumed; 


tQ  =  36.5  sec 
u  =  762  ft/sec 
v  =  0  ft/sec 
w  *  0  ft/sec 
p  -  0  rad/sec 
q  si  0  rad/sec 
r  =  0  rad/sec 
X  =  31,500  ft 


Y  »  47.78  ft 

Z  -  -13,420  ft 

t>Q  =  0  rad 

«  -G.379  rad 

^  =  0  rad 

X..  =  0  rad 
14 

=>  0  rad 


Lee  the  step  size  be  0.005  second,  the  print  interval  be  0.5  second, 
and  let  the  maximum  time  (TMAX)  be  70  seconds.  Because  the  print  inter¬ 
val  is  0.5  second  there  should  be  no  more  than  140  printouts  (MAXPT  =  140) 
Select  the  standard  printout  and  Runga-Kutta  integration.  The  input  data 
card  set  will  then  bs  as  follows: 


5  10 

I  i 

8  1 


20  30 

I  i 
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\ 

HOMING  PROGRAM 


SO 

1 


J 


70 

I 


SC 

i 


3.  Output 


The  following  is  the  output  result  from  the  subprograms  given 
in  Paragraph  1  of  Section  V  and  input  data  given  ir.  Paragraph  2  of 
Section  V. 
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Appendix  A.  QUATERNIONS 


1.  Introduction 

There  are  currently  three  methods  in  general  use  for  generating 
coordir.ate  transformations  used  in  6  degrees  of  freedom  flight  simula¬ 
tions:  Euler  angle  rotations,  direction  cosines,  and  quaternions.  This 
appendix  will  deal  with  quaternions  and  Euler  angle  rotations. 


2.  Euler  Angle  Rotations 


Euler's  theorem  states  that  any  real  rotation  may  be  expressed 
as  a  rotation  through  some  angle  about  some  fixed  axis.  That  is, 
regardless  of  what  the  rotational  history  of  the  body  is,  once  it  reaches 
some  orientation,  that  orientation  may  be  specified  in  terms  of  a  rota¬ 
tion  (through  some  angle  which  can  be  determined)  about  some  fixed  axis. 
The  theorem  can  be  restated  in  terms  of  matrices  as  follows:  for  every 

orthogonal  transformation  matrix  A,  there  exists  some  vector  R  such  that 

A  R  =  R  (A-l) 

Equation  (A-l)  is  a  statement  that  there  exists  a  vector  coincident  with 
the  axis  of  rotation  that  is  not  changed  in  magnitude  or  direction  for 
every  orthogonal  transformation  matrix  A.  That  is,  the  existence  of  an 
axis  of  rotation  can  be  established  for  any  orthogonal  transformation 
matrix  by  proving  Equation  (A-l).  The  proof  of  Equation  (A-l)  can  be 
based  on  the  more  general  equation 


A  R  =  \  R  (A-2) 

where  \  is  a  scalar  quantity  called  an  eigenvalue.  The  eigenvalue  pro¬ 
blem  is  well  known,  and  it  can  be  shown  that  the  characteristic  or 

2 

secular  equation  for  real  orthogonal  matrix  must  have  a  root  \  =  +1. 
Euler's  theorem  shows  that  it  is  possible  to  express  any  rotation  (or 
orthogonal  transformation)  as  a  single  rotation  about  some  axis  and, 
thus,  it  is  possible  to  make  use  of  the  equivalent  rotation  to  specify 
orientation. 

An  orthogonal  transformation  matrix  will  now  be  developed  using 
Euler  parameters.  Consider  Figure  A-l.  X,  Y,  and  Z  are  the  inertial 
axes  and  x,  y,  and  z  are  the  body-axis  axes  which  are  assumed  to  be 


2 

Gantroacher,  F.  R.,  Theory  of  Matrices,  Chelsea  Publishing  Company, 
I960,  Library  of  Congress  Catalog  Card  No.  59-11779. 
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fixed  to  a  body  rotating  in  space.  The  task  here  is  to  derive  the 

transformation  which  will  take  the  vector  I  =  (X,  Y.  Z)  into  the  vector 

B  =  (x,  y,  z) .  By  Euler's  theorem,  it  is  known  that  an  axis  of  rota¬ 
tion  exist  (Figure  A-l)  such  that  the  coordinate  system  (X,  Y,  Z)  can 
be  "rotated"  about  the  axis  and  be  made  to  coincide  with  the  (x,  y,  z) 
coordinate  system  for  any  possible  orientation  of  the  (x,  y,  z)  coor¬ 
dinate  system.  Now  let  che  axis  of  rotation  make  angles  a,  f3,  and  y 
with  the  X,  Y,  and  Z  axes,  respectively,  as  shown  in  Figure  A-l.  Note 
it  is  a  geometric  consequence  that  the  rotational  axis  makes  the  same 
angles  a,  P,  and  y  with  the  x,  y,  and  z  axes,  respectively.  This  fact 
can  be  appreciated  by  the  artifice  of  considering  the  X,  Y,  Z  coordinate 
system  shown  in  Figure  A-l  to  be  constructed  of  wire  and  welded  at  the 
origin  to  a  wire  representing  the  axis  of  rotation.  Clearly,  the  wire 
representing  the  axis  of  rotation  may  be  rotated  until  both  coordinate 
systems  are  coincident.  The  obvious  conclusion  is  that  tv-  es  must 
be  the  same.  The  next  step  is  to  define  a  new  coordinate  «ys .  n  (X^, 

Yr’  Zr^  suc^  xr  aligned  with  the  axis  of  rotatic  nd  such  that 

the  Yr  axis  is  in  the  X-Y  plane.  There  is  an  orthogonal  transformation 

matrix  A  such  that 


(A-3) 


In  an  identical  manner,  a  coordinate  system  (x^.,  y^,  z^,)  can  be  defined 
for  the  body-axis  coordinate  system  such  that  the  x^  axis  is  aligned 
with  the  axis  of  'otation  and  such  that  the  y^  axis  is  in  the  x-y  plane. 
Then 


* 

* 


(A-4) 


where  the  matrix  A  in  Equation  <,A-4)  is  identical  to  the  matrix  A  in 
Equation  (A-3),  Now  make  the  following  definitions: 


B  =  ( x  .  y  ,  z  | 
r  \  r  7r’  r) 


(A-4a) 
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?r  *  (V  V  Zr) 


(A-5) 


then 


B  =  (x,  y,  z) 
I  =  (X,  Y,  Z) 


(A-6) 

(A-7) 

(A-8) 

(A-9) 


Thus,  an  orthogonal  transformation  from  the  X,  Y,  Z  coordinate  system  to 
the  x,  y,  z  coordinate'  system  can  be  found  by  finding  the  matrix  A  in 
terms  of  the  angular  parameters  a>  (3,  and  7.  Because  the  X^,  axis  is 
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aligned  with  the  axis  of  rotation  and  because  the  axis  is  perpendi¬ 
cular  to  the  Z  axis  (a23  =  0) >  A  is  of  the  form 


A  = 


"cos  a 

COS  P 

a21 

a22 

a31 

a32 

cos  7 
0 


a 
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Since  A  must  be  orthogonal  matrix,  it  is  possible  to  deduce 


(A- 15) 


A  = 


cos  a  cos  p  cos  y 

±cos  p  esc  y  ±cos  a  esc  y  0 

±cos  a  cot  y  ±cos  p  cot  y  ±sin  y 


(A- 16) 


The  ambiguities  in  sign  may  be  removed  by  geometric  considerations  with 
a  =  0  which  implies  y  =  f3  =  tt/2.  The  result  is 


A  = 


cos  a  cos  p 

-cos  p  esc  y  cos  a  esc  y 


cos  y 
0 


j^-cos  a  cot  7 

-cos  p  cot  7  sin  y  J 

• 

(A-17) 

The  desired 

1  transformation  A  *  R 

A  may  now  be  computed  and 

.2  2  „2 

S  I 

2(|*|  +  £X) 

2(|£ 

-  rjX) 

A"1  R  A  = 

2(|t)  ~ 

2  2  2  2 
-i  +  n  -  r  +  * 

2(lS 

+  ?*) 

2(H  +  ■.*) 

2(t5  -  |X)  -|2 

2 

-  n 

,  ,2  ,  y2 

(A-18) 

where 

|  =  cos  a  sin  p/2 

t]  =  cos  p  sin  p/2 

£  =  cos  7  sin  p/2 

X  =  cos  p/2  . 

(A- 19) 
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The  four  parameters  g,  rj,  £,  and  X  are  known  as  the  Euler  parameters. 

It  should  be  noted  that  the  four  parameters  are  not  independent  because 

2  2  2  2 

g +  t  +  t)  +  XT  =  1  .  (A-20) 

Given  g,  £,  q,  and  X  it  is  possible  to  compute  a  unique  transformation 
matrix.  It  is  also  possible  to  solve  for  the  Euler  parameters  given 
the  transformation  matrix,  but  several  correct  solutions  may  exist, 
that  is,  the  solution  is  not  unique. 

3.  Quaternion  Parameters 

Quaternions  is  another  four- parameter  method  to  generate 
orthogonal  transformations  and  is  the  method  used  by  the  computer  pro¬ 
gram.  The  quaternion  q  is  by  definition  of  the  form 

q  =  e0  +  el  i  +  e2  j  +  e3  k  (A-21) 

where  e^,  e^,  e^,  and  are  real  numbers  and  the  vector  indices  i ,  j,  and 
k  are  defined  by 

i2  =  -i  ij  =  -ji  -  k 

j2  =  -l  jk  =  -kj  =  i 

k2  =  -1  ki  =  -ik  =  j  .  (A-22) 

The  conjugate  of  q  is  defined  to  be 

q*  =  eQ  -  iex  -  je2  -  ke3  .  (A-23) 

It  can  be  shown  from  the  definitions  above  that 

qq*  =  q*  q  =  eQ2  +  e,2  +  +  e32  .  (A-24) 

If  qq*  =  1,  then  q  is  known  as  a  versor. 

The  quantity  eQ  is  called  the  real  or  scalar  part  and  ie^  +  je„  + 
ke3  is  called  the  complex  or  vector  part.  Now  let  V  be  a  quaternion 
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whose  scalar  part  is  zero.  Thus  V  may  be  thought  of  as  a  vector 


V  =  iX  +  jY  +  kZ 


(A-25) 


Now  consider 


q*Vq=V' 


-  (A-26 ) 


where  q  is  a  versor  (q*  q  =  qq*  =  1). 

It  can  be  shown  that  V'  is  a  , , 
lished  for  quaternions,  and 


by  using  the  definition  estab- 


V'  =  (eQ  -  ie.  -  je2  -  ke^  |iX  +  jY  +  kZ^e0  +  iex  +  je2  +  ke^(A-27) 
expanding 

V  =  x|X  [sq  +  e^  e2  e3  ]  f  Y[2e3e0  +  2ele2  J  +  Z  [2ele3  “  2e0e2]} 

t  jjx  [2eie2  “  2e3eo]  +  Y[e0  "  ei  +  e9  "  e3  ]  +  Z  [2eleC  +  2e3e2]| 

+  k{XK62  +  2elC3]  *  Y[2ele3  '  2Vl]  +  Z[£02  “  el2  "  e22  +  e32]} 

(A-28) 


or 


V'  = 


2,  2  2  2 

e0  **1  "e2  "e3 


2(elVe3eo) 

^‘iWz) 


'(e3e0+ele2/ 


2  2  2  2 
e0  "el  ~e3 


-(ele3'e0e2) 

!(ele0+e3e2) 


2(e2e3_e0el)  e02"el2-e22^3 


(A-29) 


The  aoove  matrix  can  be  shown  to  be  ar-  orthogonal  transformation  matrix. 
Note  the  3X3  transformation  matrix  ie  completely  defined  by  the 
quaternion  q. 
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The  above  matrix  is  used  to  transform  vectors  from  the  inertial 
system  to  the  missile  body-axis  system.  There  is  a  one-to-one  corres¬ 
pondence  between  the  matrix  in  Equation  (A-18)  and  the  quaternion  q. 
The  following  relationships  may  be  established  by  comparing  Equations 
(A-18)  and  (A-29) 


e 


0 


X 


-  n 


This  correspondence  shows  that  the  set  of  all  matrices  of  the  form 
of  Equation  (A-29)  is  the  same  as  the  set  of  all  transformation  matrices 
obtained  by  Euler  angle  rotations. 

The  standard  Euler  angle  ft,  $  (Figure  1)  can  be  computed  from 
the  quaternions  by  the  following  relationships: 


f2  (*1*3  -  V2)) 
/  2(el62  *  V3)  \ 

W  +  *i2)  -  v 

/  2(e2e3  *  Vl)  \ 

W  *  V)  -  »/ 


(A-30) 

(A-31) 

(A-32) 


which  can  be  established  by  comparing  the  transformation  matrix  using 
quaternions  to  the  well  known  three  parameter  Euler  angle  transforma¬ 
tion  matrix.  Similarly, 


=  cos  Ifi 2  cos  6/2  cos  if 2  +  sin  v/2  sin  0/2  sin 

=  cos  y/2  cos  0/2  sin  tf>/2  -  sin  '1/2.  sin  9/2  cos  $/2 

=  cos  v/2  sin  ft/2  cos  i  f 2  +  sin  %/2  cos  ft/Z  sin  */2 

*  -cos  v/2  sin  0/2  3in  $/2  +  sin  if 2  cos  ft/2  cos  b/2  .  (A-33) 


4 .  Cayley-Klsin  Parameters 


This  section  will  lay  the  groundwork  for  the  next  section 
where  the*  relationship  between  the  angular  rates  p,  a;  and  r  r.nd  the 


quaternions  will  be  developed.  The  basic  idea  of  the  Cayley-Klein 
parameters  is  to  represent  a  real  rotation  (or  transformation)  with  a 
2X2  complex  matrix  instead  of  the  usual  3X3  real  matrix.  Thu3.  analytic 
operations  are  great?y  simplified.  Consider  the  following  complex 
matrix 


H  = 


lll 

h22 

21 

h22 

.< 

(A-34) 


Let  the  matrix  have  the  following  properties: 

a)  H  is  a  unitary  matrix 

b)  |  H !  =  -fi. 

Then  from  tho  unitary  property 


H*  H  «  1  »  H  H* 


Thus 


k 

hi 

*  1 
h21 

1 

1 

* 

r 

[*u 

h12 

1 

0 

■>7 

h12 

* 

h22 

I 

h21 

^22 

0 

m 

1 

Expanding  and  equating  coroosiants 

hil  hll  +  “21  n21 


=  1 


hn  h12  *  h21  h22  =  0 


h12  hll  +  h22  h21  =  0 


hlZ  h12  h22  h22  *  1 


(A-35) 


(A- 36) 


(A-37) 

(A-38) 

(A-39) 

(V40> 


4S 


Note  that  the  left  hand  side  of  Equations  (A-38)  and  (A-39)  are  the 
complex  conjugates  of  each  other,  thus,  there  are  only  3  Independent 
equations.  Now  from  the  fact  that  | H j  =  +1,  we  have 


nll  h22  “  h21  h12  "  41 


(A-41) 


From  Equation  (A-38),  we  have 


(A-42) 


I 


Then  using  Equations  (A-41)  and  (A-42) 


1 


(A-43) 


Since 


h0^)  =  1  from  Equation  (A-37) 


(A-44) 


Similarly 


h 


22 


(A-45) 


So,  the  matrix  K  may  be  written 

I”  hll  h12  ] 

(A-46) 

The  quantities  h^,  h^t  h^1 .  and  h,,9  are  usually  referred  to  as  Cayley- 

Kiein  parameters  and  are  in  general  complex  numbers  which  con  be  defined 
by  the  two  complex  numbers 


H  *  1 


I  "hi2 


11 


l 

J 


4 

1 
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•ilAWAx kw*#* 


hll  =  co  +  lci 


®  ^3 


and  the  matrix  H  may  be  written 


+  iCj 

C2  + 

iC. 

+  iC3 

C0  " 

iC. 

Now  consider  a  matrix  P  which  has  the  following  form 


x  +  iy 


x  -  iy 
-z 


where  x,  y,  ^  are  real  numbers  which  can  be  viewed  as  coordinates  or 
components  of  a  three-dimenc.  ion-1  vector  (X,  Y,  Z)  in  Euclidean  space. 
Since  P  is  Kerm'tian,  the  transpose  of  the  complex  conjugate  of  P  is 
equal  to  the  matrix  P  and 


P*  -  P 


Now,  consider  the  similarity  transformation  of  P  of  the  form 


P:  =  H  P  H*1 


The  following  properties  of  P  are  iovarient  under  a  similarity  trans¬ 
formation:  Kermitian  property,  trace,  and  determinant .  It  follows  P 
must  have  the  following  form 


it  follows 


2 

x 


+  z 


x'2  +  y'2 


+  z 


,2 


(A-54) 


If  x,  y,  and  z  are  components  of  a  vector,  then  the  length  of  the  vettbr 
has  not  been  changed  by  the  similarity  transformation.  Then  from 
Equations  'A-49),  (A-51),  and  (A-48) 


z'  x'  -  iy'" 

‘C0+1Ci  C2  +  1C3‘ 

N 

M 

1 

C0  -  icx  -C2  -  ICj 

x*  +  iy’ 

m 

*c2  +  ic3  CQ  -  iC, 

x  4  iy  -x 

S  *  “3  c0  +  IC1 

m  a 

. 

U 

(A-55) 


Then  from  Equation  (A-55),  it  can  be  shown 


”  x'“ 

'4  -  cj  -  c\  +  4  2(coCl  +  C2C3)  2(ClC3  -  C0C2) 

y' 

m 

2(C2C,  -  C0Cj)  Cl-Cl  +  C\  -  4  2(ClC2  +  C3C0) 

2(C0C2  +  C1C3>  2(C1C2  -  CoH)  C0  -  Cl  -  C2  '  SJ 

"  x  ■ 

y 

z 

(A-56) 


Using  matrix  notation 


(A-57) 


The  matrix  A  satisfies  orthgonality  conditions.  The  parameters  Cq,  C^, 
C2*  and  may  be  related  to  results  in  the  previous  two  sections  by 
comparing  Equations  (A-56),  (A-29),  and  (A-18). 


(A-58) 
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Thus,  an  equivalence  ha,  been  indicated  between  the  real  3X3  matrix  A 
and  the  complex  2X2  matrix  H. 

»  -  ~  zt 
J£.“£EEV£  ««»  w  «■ 3X3  ^ 

*•  .  B  *  ‘A-59> 


Now  let  the 


associated  2X2  complex  matrix  be  Hj,  so  that 


p^^rh'/ 


(A-60) 


Consider  a  second  transformation  A  associated  with  H2 


r"  =  A  r 


Now  substitute 


P "  -  S2  P  H? 

Equations  (A-59)  and  (A-60)  into  Equation  (A- 61) 


r"  a  A  B  r 

?' '  =  H2  Hj  p  k"1  h'1 


(A-61) 


(A-62) 

<A-63) 


but  A  B  .  C  and  it,  H,  -  H.  uhere  C  ia  a  real  3X3  mtrlx  and  H3  is  a 

ccplax  2x2”Jtl'fLh”^-.l'3X3f0»™riLs‘l''or«sP™dfto  ttaT^UpUcation 
3lS1SSSt 2  complnxhacrices  in  t«  sa«  order.  Thus  there 

“r/reirarael^i  rratrix. 

£„rMfi™r»atrL?iih”bruLf  t‘f  dehv;5Sr“?respondancebet»eenthe 

Aguiar  rates  p.  q,  and  r  and  the  quaternions  V  V  «r  «3 
next  paragraph. 
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The  transformation  matrix  using  quaternions  has  been  developed 
in  the  preceding  section.  Since  the  body-axis  angular  rates  p,  q,  and 
r  determine  the  orientation  of  the  body-axis  system  relative  to  the 
inertial  system,  the  relation  between  the  rate  of  change  of  the  quater¬ 
nions  (4q,  e^,  c,,  e^)  and  the  angular  rates  (p,  q,  and  r)  must  be  estab¬ 
lished.  It  was  shown  in  paragraph  4  that  an  orthogonal  transformation 
matrix  may  be  represented  using  the  Cay ley- Klein  approach  by  a  2X2  matrix 


C0  +  iCl 


C2  +  iC3 


-Cg  +  iC^  Cq  ”  iC^  • 


(A-64) 


Then  using  Equations  (A-64)  and  (A-58) 


cos  -r  +  i  cos  7  sin  ^ 


cos  p  sin  ^  +  i  cos  a  sin  ^ 


-cos  p  sin  —  +  i  cos  a  sin  ^ 


cos  ^  -  i  cos  y  sin  ^ 


(A-65) 


Now  let  u  =  A  u  be  an  infinitesmal  rotation.  Then  if  we  assume  cos 
Au/2  =  land  sin  Au/2  =  au/2, 


i  ,  Au 
1  +  -p  cos  7 


A  n  .  A  n 

— p  COS  P  +  1  -p  COS  a 


=p  cos  p  +  i  ^  cos  a 


.  An 

l  -p  cos  7 


(A-66; 


Now  assume  that  the  rotation  Au  occurs  during  the  time  At.  If  H  is  the 
matrix  at  the  beginning  of  the  rotation,  then  H  is  the  matrix  at  the 

end  of  the  rotation,  and  the  time  derivative  of  H  may  be  written 


then 


i  cos  7 

dH  _  1  dji 
dt  2  dt 

_-cos  p  +  i  cos  a 


Because  d^/dt  is  the  scalar  quantity  of  the  angular  velocity  vector, 
the  p,  q,  and  r  components  of  this  velocity  vector,  as  defined  in 
Figure  3,  are  given  by 

dp 

p  =  dt  COS  a 

dp 

q  =  dt  cos  p 
du 

r  =  —  cos  7  .  (A-69) 


2  C0  =  C3P  '  C2q  '  V 
2  C1  =  -C?P  +  C3q  +  CQr 

2  C2  =  ClP  +  CQq  -  C3r 

2  c3  =  CQp  -  C^q  +  C2r 


(A- 71) 
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Then  using  Equation  (A-58),  the  relationship  between  quaternions  and 
angular  rates  is 

40  -  -1/2  (ejP  +  e2q  +  ef) 

H  *  1/2  (V  '  e3q  +  e2r) 

42  -  1/2  (V  +  e0q  -  V) 

e3  =  1/2  ^-e2p  +  e3q  +  eQr\  .  (A-72) 

It  will  be  recalled  from  Paragraph  3  that  tha  quaternion  q  =  e^  +  ie^  + 
je2  +  ke^  is  a  versor  which  implies 

e0  +  el  +  e2  +  e3  =  1  ’  (A- 73) 

Equation  (A-72)  is  used  in  the  computer  program  to  compute  the  rate  of 
change  of  the  quaternion  components  subject  to  the  constraint  given  in 
Equation  (A-73).  Hie  matter  of  a  constraint  is  handled  computationally 
by  rewriting  Equation  (A-72)  in  the  fora 

eQ  =  -1/2  /e,p  e2q  +  e^  +  K  e 

41  *  1/2  (°0P  -  e3q  +  V)  +  K  *  81 

e2  =  1/2  ^e3p  +  eQq  -  +  K  e  e2 

e3  =  1/2  (~e2p  +  e.q  +  e^r^  -f  K  ?  e3  (A- 74) 

where 

e  =  1  -  (ej  +  ej  +  e \  +  e*}  (A-75) 

and  where  K  is  an  arbitrary  real,  positive  constant.  The  program  uses  a 
value  of  K  =  100  which  was  found  by  empirical  methods  to  be  satisfactory 
for  all  cases  tested.  The  value  of  K  may  be  modified  in  the  program  by 

reading  it  in  through  one  of  the  subroutines  supplied  by  the  user.  How¬ 

ever,  large  values  of  K  should  be  avoided  because  the  result  will  be  an 
unstable  solution. 


Appendix  B.  INTEGRATION  ROUTINES 


1.  Introduction 

There  are  three  integration  routines  supplied  with  the  main 
program:  Runga-Kutta,  Runga-Kutta-Merson,  and  Hamming  Predictor  Cor¬ 
rector.  The  user  selects  the  desired  integration  routine  by  an  input 
data  card.  The  step  size  for  all  three  routines  may  be  altered  by  the 
user  at  any  point  during  the  simulation  by  redefining  the  FORTRAN  vari¬ 
able  DT  in  a  user  supplied  subroutine  which  contains  the  common  block 
0THER2.  It  should  be  noted  that  all  the  integration  subroutines  have 
identical  arguments  in  the  call  statements  but  that  some  of  the  FORTRAN 
variables  in  the  argument  lit  t  may  be  dummy  variables .  The  equations 
to  be  integrated  are  generated  by  the  EXTERNAL,  SUBROUTINE  DESUB  (.TIME, 
X,  DX). 


2.  Runga-Kutta  Integration 

The  Runga-Kutta  integration  routine  is  a  fourth  order  method. 
The  call  statement  is  of  the  form 

CALL  RUNGA  (TIME,  X,  DX,  R,  DT,  DTMIN,  EMAX,  NT,  IC ,  SMAX,  DESUB)  . 
The  set  of  equations  to  be  solved  is  of  the  form 

xi  -  fl(Xl'  X2 . V  c) 

X2  *  f 2  ’  X2 . V  C) 


Xn  *  £n(Xl’  X2 . V  ')  '  (B-1) 

Then  the  equations  used  in  fixed  step  Runga-Kutta  integration  are 

X(tK  +  l)i  -  X('k)i  +  1/6  (Kil  +  2  K12  +  2  *t3  +  Km) 

1  *  1  j  2  j  *  ♦  •  ^  n 
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where 


a  - h  •  fi(xi(sc)-  -  •  •  •  x„(£k)-  ck) 

i2  “  h  '  fi(Xl(tK)+  KU/2'  ^(V)  +  K21/2’  •  •  •  >  xn(CK 
\ 

:)  +  k»i/2;  £k 

+  h/2) 

(B-2) 

13  =  h  fi(Xl(tK)  +  K12/2’  X2^k)  +  K22/2’  *  *  *  ’ 

k)+  Kn2/2;  fcK 

+  H/2) 

14  =  n  fi(Xl(tK)  +  K13’  X2 (Ck)  +  K23’  *  *  *  *  Xn(fcK 

)  +  Kn35  lK  + 

where  h  is  the  step  size  and  t  denotes  a  specific  time  at  step  size  h 

K 


intervals  (tv  =  K.h  where  K  =  1,  2,  .  .  .)• 

K. 

This  program  has  the  advantage  that  it  requires  less  com¬ 
putation  than  Runga-Kutta-Merson  and  is  faster  when  there  are  no 
fast  transients  (relative  to  the  other  transients)  which  die  out  and 
force  the  use  of  a  small  step  size.  The  program  also  has  the  advan¬ 
tage  in  that  the  user  can  determine  what  stage  the  computation  has 
progressed  to  in  the  subroutine  by  testing  the  variable  R(4) .  The 
variable  R(4)  has  the  following  significances:  R(4)  =  1.1  implies 
K.^  is  being  computed  (starting  with  subroutine  SEEKER),  R(4)  =2.1 

implies  is  being  computed,  R(4)  =  3.1  implies  is  being  computed, 

and  R(4)  =  4.1  implies  K..^  is  being  computed. 


3.  Runga-Kutta-Merson 


The  Runga-Kutcd-Merson  method  is  a  variable  step  size  program 
which  in  certain  situations  is  potentially  faster  and  more  accurate 

3 

than  the  other  methods.  This  method  is  a  fourth  order  method  which 
uses  the  following  equations: 


martens,  H.  R.  A  Comparative  Study  of  Digital  Integration  Methods, 
Simulation,  February  1969. 
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Kil  *  1/3  h  fi  [X1^k)1  X2  (Sc)  ’  *  '  •  *  xn  (Ck)  *»  tJ 

Ki2  «  1/3  h  f±  [Xi^K)+  Ku>  X2(tK)+  K21»  *  *  *  »  Xn(tK)+  KnP  *11  +  h/3] 

Xi3  "  1/3  h  fi  [Xl(‘K)  +  *5  Kll»  +  *5  K12  >  hfy  -5  K21 

+  .5  K22,  .  .  Xn  (t^+  .5  Knl  +  .5  Kn2;  tR  +  h/3] 

X4i  =  1/3  h  ^  [Xi  (fcK)  +  3/8  K11  +  9/8  K13’  X2(tK)+  3/8  K21 

+  9/8  K23,  ....  Xn(t^+  3/8  Kpl  +  9/8  K^;  tR  +  h/2] 

X5i  "  l'3  h  fi  [Xl(tK)+  3/2  K11  "  9'2  K13  +  6  K14‘  ^(‘k)  +  3/2  K21 

-  9/2  +  6  K24>  ....  Xn(tR)+  3/2  Knl  -  9/2 

+  6  Kn4;  fcK  +  h]  CB-3) 

then 

X£(k  +  1)  »  X£(k)  +  .5  ^KiL  +  4  K.4  +  X15^ 

where  i  *»  1 ,  2 ,  .  .  . ,  n 

The  estimate  of  the  truncation  error  is 

-i  =(Kil  *  9/2  Ki3  +  Ki4  ~  1/2  Ki5)  /5  * 

The  step  size  DT  is  changed  by  the  program  until  SMAX,  the  maximum 

element  of  the  set  { |  ej  ,  1  e2l  ,  .  .  .,  je^!},  lftss  than  EMAX,  subject 

to  the  restriction  that  DTMN  5  DT  where  DTMIN  is  the  minimum  allowable 
step  size.  If  the  computed  DT  is  less  than  DTMIN,  DT  is  set  equal  to 
DTMIN.  If  SMAX  S  EMAX  for  three  steps  in  a  row  the  step  size  DT  is 
doubled,  and  if  SMAX  ^  EMAX  the  step  size  DT  is  cut  in  half  (DT^  DTMIN). 

The  calling  statement  for  the  Runga-Kutta-Merson  integration 
routine  is  of  the  form 

CALL  RKMER(TIME,  X,  DX,  R,  DT,  DTMIN,  EMAX,  NT,  IC,  SMAX,  DESUB). 

The  parameter  R(4)  has  a  definition  which  is  very  similar  tc 
the  parameter  R(4)  in  the  Runga-Kutta  routine.  The  variable  R(4) 
has  the  following  meaning:  R(4)  «  1.1  implies  X£l  is  being  computed. 
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R(4)  =  2.1  Implies  is  being  computed,  R(4)  =3.1  implies  K^3  is 
being  computed,  R(4)  =  4.1  implies  is  being  computed,  and  R(4)  s 
5.1  implies  K^,.  is  being  computed. 

4.  Hamming  Predictor  Corrector  Intergration 


This  integration  routine  is  a  modified  version  of  IBM  SUB- 
4 

ROUTINE  HPCG.  The  subroutine  is  basically  the  same  except  that  it  was 
changed  to  be  compatible  with  the  other  integration  routines. 

Hamming’s  modified  predictor  corrector  method  is  a  fourth  order 
method  using  4  preceding  points  for  computation  of  a  new  vector  of  the 
dependent  variable  X.  A  fourth  order  Runga-Kutta  method  is  used  for 
adjustment  of  the  initial  step  size  and  for  computation  of  starting 
values.  The  routine  automatically  adjusts  the  step  size  DT, halving  or 
doubling.  If  the  step  size  DT  is  halved  more  than  10  successive  times, 
an  error  message  is  generated. 

The  calling  statement  is  of  the  form 
CALL  HAMPC  (TIME,  X,  DX,  R,  DT,  DTMIN,  EMAX,  NT,  IC,  SMAX,  DESUB) . 


^IBM.  System/360  Scientific  Subroutine  Package  (360A-CM-03X) 
Version  III  Programmer’s  Manual,  H20-0205-3,  1968. 
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Appendix  C.  MAIN  PROGRAM  LiSTING 


C  THIS  °ROGRAH  IS  A  GENERAL  PJRPGSt  PROGRAM  FOR  THE  6C  SIMULATION  OF 
C  TERMINAL  HOMING  MISSILES 

C 

G  COORDINATE  TRANSFORMATIONS  ARE  MADE  WITH  QUATERNIONS 

C 

C  —  THE-  USER  -MUST  5UPPLY  -THE  FftfctGWI-NG  SUBROUTINES-  - -  - -  -  ■ 

C  ScF.KER 

C  TARS ST  -  •  - 

C  FORD M 

0  -MRT  -  -  •  ... 

C  STRFLf 

C 

&  USER  INSTRUCTIONS  FOR  THIS  PROGRAM  MAY  3E  FOONC  IN  R£-TR-72-i6 
C  21  SEPTEMBER  1972*  US  ARMY  MISSILE  COMMAND,  REOSTONE  ARSENAL. BY 

€-  OR.  LEWIS  G.  MINOR 

C  OR ,  LEWIS  G.  MINOR 

C 

je*», 

C  *  LIST  OF  VARIABLES  USED  IN  PROGRAM  * 

£*****{■»*  *.  ««  J.»  «  *  -i.  .»  4  .*  It  *  f  *  *  «  t  »  W  •  *.  «.  M  .*  .X 

c  x;d  =  p 

e--x-'2-)»o— . . . —  -  - - 

C  X{3)=R 

C  X<<rT*£-3 .  -  -------- 

C  X  (55  =£1 

C  X(0)=E2  -  -  -  - 

C  x  (7) =£3 

6  -X{flf-*F.-&9T - - 

C  XO)=Yfc  COT 
C  X(lj,=7-  hot 
C  X(11U,!E 

£ Y  —  yr .  ..  _ _ _ _  ..  .  __  .  _  .  _  _ _ 

x  x  c  t  *  .  —  *  • 

C  X ( 1 3) ~ZS 

G  «MI  -  1-1  4,-5  M-ii-  —  ARE — THE  -SEE*  £-R  -  STATES -  - . -  - 

C  X(J>  J=14>SX,SXfT**i3  ARE  THE  TARGET  STATES 

0  -  ....  . 

C 

DIMENSION  A(53>TOXf£W  --  -  -  -  -  — 

OIMENSION  TAGC  8> 

EXTERNAL  .S££y*~  - -  -  - -  - . — 

REAL  MACH 

REAL-  -MXtMY-,MZ-tK,K£  — .  -  -  •  -  - 

REAL  MASS 

REAL  IX.IT-vJ  ?,  IZX  .  _  .  - 

INTEGER  SX.TX 

■COMMON  -/TRANS/ - -  -  -  --  —  -  - 

COMMON  /OISPL/  U,V,W,OELXE,p£LYE,DELZE,XTE,YTE,ZTE,RM 
COMMON  /ANG/-AL-FA-,8ETA-,-SIGY,5  I-GZ,  0R1,DY2,QP3,DY4 
COHHON  /AIR/  HACH.VH, VS,QS,QSO ,QSX,RHO 
COMMON  /FRHO  /  FX ,  FT, FZ, MX,  MY , MZ , IX , 1 7, 1 Y , IZX ,NAS3 
COMMON  /OTHERl /ITERA .DT.OTMIN, EM A  X.SMAX.TMAX, PRNTI 
OONMHN/QT+ER  2/  TH?ATR-SI-,RHlTKE-rS  ,-0t3X-,TX-,K,G,  KAXRT 
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-  4 I-HJ/  Sri*  !- - 

CONHON/MIS/  RMIN 

1-  -  IV£R4=i  -  . . 

RMIN -109 J. 

TRF=9 

lKftiCPT=0 

IPLOT=T 

--  -IXI-t  . 

C 

c  _*  -5  .*  3L  4  -5L  Jf_4.  i  Jf..*  .*_*  x  1  *-J  ..»  .«■.**  »  ♦  *  .»  * 

C  *  READ  INITIAL  CONDITIONS  * 

M- - jrC-S^yjjCT^  ^Ut)OAa&-m«WU  lg,.lQgTM*tT~S»fcftXAk-ORWWA - * 

C*  ir  ICi?TN=2,STAN0AR0*-SP£CIAL  PRINTOUT/  * 

- - - 

O' IF  INTOPT=C,  FIXED  STEP  RUNGA-XJTfA  INTEGRATION  USEO  * 

-  ■  ~iF  ^Ton.T-=4^4^Nefc^CT-TA-NE^^aU-  I  MI  £RSRATXQM--*FSEU - * 

C*  !T  I NTOP  f =2*  HA  THING  PREDICT  OR  CORRECTOR  INTEGRATION  USES  * 

-G-2 -  *-.*_*_  9  9*9-  £_£  _t.JL  Jt_»_ JL..TL.JL-T— -t—TL-f  JUS.  .*....? 
RSA0(5, 201}  SX»TX 

..  .  REAO-CS,5WJ-OT-  it  ASTI) ,  T=W&1  -  _ .  .  . 

800  FORVATISAIG) 

XFISX.C^.tT  -CALL  EXIT. .  ..  .. 

READ  15,201}  IOPTN,INTOPT 

NSXl-NSX+l 

NTSX-T'CtSX  .  .  -  .  —  .  . 

NT=SXFTXT13 

REAOiStii*  TIME,  PHI,  THTA^PS-t  ... 

READ  15,11)  U,V,R,XI1) ,XI2>,/(3) 

...  —  «EAS15^VT.Xtiil.-XU2},X-U31 . . -  - - - 

REACT?, US  IX  Til 

REACT?, Ul-  -  U  UUUUSXt'Kti .  .... 

C 

£*»».***•**,>*****»■**•»***■*<»**#*  *.««■»» 

C  «  READ  STEF  SIZE,  PRINT  INTERVAL,  AND  HAX  TINE  • 

c  -T  ?-  V  4.  *  ’.. _1.  7.  JL  *_*_  *  »  »  >__*  . »  .».».♦•»»»»»» 

RE 43  *5,11)  OT,OTHH{,EHAX,PR'»TI,THAX 

RLA!)  <.5,2011.  HAXPT  .  ....  - -  .  — -  - 

c 

Q  »  *  *  .?  ■»  ,  *  *  *  *  ■•  *»  «  c  *.*.**  3  ******  «■***_** 

C  *  OTHER  CONSTANTS  * 

C-  V  »  f.  ♦  *»»»«  *  *_*  *.*_.*_?_*  *..*  *  *  *  *_  *  *  *_*  *  J»  *  *  .f  *  ♦  *, 

C-32  «  1 1'49  A 

XJUTPTrrlthE-  .  ..  -  .  .  - . . 

K=10  9  . 

C 

£.*.***,*  _*■_*_  V  *  *  *  *  T»r«  ********  *  *  f  *  *  •*_ 

C  *  COMPOTE  IC  rOR  E& ,Ei .fcZ, E J,0XE ,OYc ,OZE  * 

C  *  4  *  *  *_  *  >•  f  *  i.  *_  *  *  v  *  *_*_  .*„_*  *  *  *  *  _*  ♦  »  e_  _*  ^  _o_  *  _*_  * ;_X 

C0SPSI=C0STPS£/2U 

.  -U0STaT=C0SUiiTJT/2^»  -  - -  -  - - - - 

COS?HI=COSTPHI/2.) 

- -smF-s!=s-i*up-sixu x - - —  —  ...  - - 
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C  EO 


—  3-fNT  H  T-*5 1 U  vERF-A/gr* — 
SI  N?  MI-SIMP  HI /2.  > 


X  (4)  =C'Oi.‘P?,I',CJSrHr*COSHH;<-SINJ;>S:*SINTHT*SINPHI 

€  Cl  . -  -  -  -  - 

X  (5;  =C0SPSI*C3STHT*SINPHI-SINPSI*SINTH?*C0SPHI 

X?&) =COSPSI*SI NTKT*C  CSPPI+S i NPSI *C03THT*SINPHI 

C  £3  .  - 

K  17) =-CD3PSI*SINTHT*SIKPHI+SINPSI*C0STHT*C0SPHI 
GH.L-  feFT0E*f*#-X<4),X<5>,Xl6),X{7)>  - 

CALL  MULT (X<  8?  ,X { 9)  ,  X i 10 ) ,ArU , V, K) 

- 6-ftfct— W: - 

R(i)  -QT 

-G0  TO -51-3  -  -  -  - 

C 


€ 

C 

9 


»-  9—  *  *  »-*  -  *  9_  *  *:  *  *. 

*  start  integration  l*jcp 


*_•»  *  »  « 


CONTINUE 

ITERA^O-RA-H 


9  *  .  J» 


*  *  -  9  ■*  *  *.  *. 

-  9.  _*  -C  -  » — • — S. 


If  aNfOPT.Ca.l  J  GO  TO  30 
IFdFM'OPT.eQ.Z)  GO  TO  31 

CALL  RUNGMTIfiE.r  ,PK  ,R,0T,0rMll'.',EMAXtNT,IC,SMAX,0£SU8) 

3 5/  CALL  (?KMER'riMEtX,DX,R,OT»CiTMIN>  EMAX,NT,IC,SMAX,0ESU8> 

--  GO -TO  32-  -----  -  .  -  -  - 


31  CA»u  HAMPCiTIME,X,0X,R,0TfDTHINrEMAX,Nr,IC.SHAX,DESU3) 

32  CONTINUE-  -  •  -  - 

£**»*******«**♦♦<>*  ?**?**«*****•»***•? 

-G  * Pft©Of^A^-PR4N-f-OU4 - — -  - - 

IFITINf.LT. XASTPT.ANO./TERA.NE.D  GO  TO  41  - 
513  XA^T°T=/AST°T>  PP.H  T I 


-  -  -  IM*OP-TW,£a^»— £<J  TR-&>*- . . - 

IF (I TER6  »  ME. 1)  GO  TO  512 


. -MKS T-g-<frr-»R34— <  TAG <  I>-T4>1 .  S) - - - - 

831  FORMAT (1H1.H0X .  3A10X//) 

WR-fF F<6,2t4J UT-,3?MIN-TE^A^.TNAX-rPRNTI-  ■  - . 

512  WRITE  !6, 13)  TIME 

HPIT€<6,144  tfMrKACMTXM13,<(12F7XF13> 

THT4  =  4SIN  (  -2  •*  (X  <5}*X  (75  -X<4>  +Yib))) 

-  -P-51-4  TAN2-f3»-*F-X-‘-5  )  *X  < 6T  ■> M7T-)-r^x-MX'-4-MX-F-41-S  X45T*X  <5U-FV) 
Pt':t=ATAN2C2.*(X{6J*xrf  >*X14I*X15)  )  .2.*  «  SM  *x  (4)  -X  {TJ  *  XT  7) ) -i  .  i 

WRITE <5,555*  -T^TArPSItP+fl .  -  -  - 

WRITE  <6->S0C)  J  ,  V»W,ALrA,9ETA 

WRIT-?  <6,  164-  3«itX€,U&LAf£rO€i.ZE.F^MX  -  -  - - 

WRITE '6.5015  X 11) ,X12>,X?3>,FT,HT 

— HRI-T  €  <  Gr»^2*-&-PyT  UT2  t-PHtF-Zt 4 i-  — . . . . . 

WRITE  10,503)  <3P3*OY4,SIGY»SIGZt<S 

IFdOPTN.FQ.  03 — -GO  TO  701  -  ...... 

633  CALL  WRT ITIMF,\.CX> 

7*3 1  IKAXPT=IP.AXPT4-1  -  -  ... 

IF<THAXPT  .LE.MAXPT)  GO  TO  41 

WRITE  •  --.‘+AX-PT-  - -  -  -  - -  --  - 
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-GG--IO  -4- —  - - 

*>1  CONI  I  \'CZ 

GALU  STRPLTT-TI-HEr-XrOX-rlMlSSf  IP40T)  . . - . 

IF(I  PLOT  t  £Q«  0)  GO  TO  637 
IFdHISSi  60  5» -605*  60  6  ■ 

637  CONTINUE 

•  -  -if -41  !*£-,  &T*T  WVX4  GO- -TO -4 - - - -  - 

C*  «*»**«* 

C*  — CGMPUTE-NJ-SS  OI-STANOE  -  -  - -  -•  -- . .  * 

r  t  *  *  *  *********  *********.*♦****** 


- Gf-4X414i-..fe6.0*4— -&0-TO  642 - - 

IFdXI.gj.e.f  GO  TO  612 

- - — *3) - 

WRITE  (6,669) 

- K-ftTT£<W6-U4--ROH.rX«t4-)-fXt.G01.f£-T430 - 

WRITE<6,9Q9) 


-549  - 
-6-44-  - 
639 


FORX-4  14 / / *X_X~X -X-X  -X-X-X-X-X-X-. X  -X  -X  X--X- -X-X  -X-X-  X -X-X-X  X  -  X-X-X-- 
1XXXXXXXXXXXXXXXXXXXXXXXXXXXXXX  x*//> 


-FORHA 14*-  4H-5T-ANGE  -FROM-  TAR-6E-T—  A-T  - 


1*  ‘  X£=*,£ii.4,5X,*YE  =*,£11.«*,5X  ,*ZE  ***E11.4//> 

FOfi))A-T(///*  X  xxx-xrxx  XX  X-XXX-XXXXXXXXX.XXX 
1EARTN  IMPACT  POINT  XX).  XXXXXXXXXXXXXXXXX  X*//5 


612 


631 


635 


1X1=0- 

IF (RMfGT  . 50. )  GO  TT 
.  i.FQ-T.GT^-WOJG)  4U-  v  ‘ 

if(rm.gt.pmin>  iff -i  ft 

IFdFF.GT.20)  C-0  TO  601 
Tr  IRWIN, i_F..RW>  GO  TO  502 
-  OXMIN*GELX£ 

0YM£n=0ELYE 
-4)Z«I*=0£42£~  -  - 

RMIN=RH 


INTOPT*a 


GO?0  602  -  -  - 

WRIT  £ (6, 302) 

WRIrE  J6.3C31  PHI M*OXHUi«QT.'iI-tt.3 ZNa.N 


WRIT  £ 36*91!!) 

-JU>VU  =  1- . . . —  ---  —  —  ..  - - 

IKISS=3 

-GO  TO  41  . .  ..  . 


CONTINUE 


GO  TO  1 

630  IF(X i 13) •  CT. 5  J •»  GO  TO  602 

X?*LOT=;  .  -  -  -  - 

XMISS=i 
-WRITE46*  3U2) 

WRITE (fi»  6G4) 

-WRITF  (6*916) -  - 

GO  TO  41 

Si)  6  -CONT-I  Wtlc  -  - - - 


GO  TO  1 

602  CONTINUE  -  -  '  -  -  - — . .  . . . 

GO  TO  9 

60J'  ?JRW4Tty/V4X*ii*rfHJSS  Ol£TANG£--GR£.aTER-DUN-£S-F'£^3  - - -  -  _) _ 

332  FORMAT(///lX  »123H-t  *•  ♦  «■  *  *  *  *  *  *  *  *  *++*♦*+*++.*.  +  + 
...  .  4 niter  0IS-7ANC 6* 4-  4 .A-A-A-A.  A -A.  A-A-A-+-A-4  .  »  A —  
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303  FORMAT  (//IX,  32  HCLOSEST  APPROACH  TO  THC  TARGET  =tEii.4,5H  FEET/ 


SUBR-OWfN€  BfWEA  <4,£«>,E4r»2-2-r€3-i 
9 1  ME  N  S  ”  .'  N  Ail) 

£.j$=i'  ■«**  J 

£1S=  E  1*E1 
E2S=E2*E2 
E3S=F3*E3  ■ 

-€i2—£  i*-E-e -  -  - -  - 

E03=ED*E3 
£13=  E 1*E 3 
Ej2=E0*E2 
E23=E2*E3 
E01=EG*E1 

A44-)  =&»S*E1>-**E  2-S^£-3S —  - - 

A  { 2)  =2.*  I  E 1 2  +E  0 3 ) 

At3^=3.M-El^E-421-  -  -  . . 

A  <4»  =2.*  (E12-SD3) 

A  <5)  =E«5  +  E2S'E  1S-E3S 
0  lb)  =  2.MF23*Ell) 

A  (7)  =?.<*  (E  J3*E  *>2) 

A  ( 8)  =  2.*  (E23-E  Cl) 

A(q)  =ECS  +  F3S-E1S-E2S 

RETURN 

END 


SUBROUTINE  EATG3FCAI,EotEl,E2»E3) 
DIMENSION  t[m,A«9) 

CALL  RBT  GE  A(  A,  E3»E1»E2»E3) 

Alii! =6(1) 

AI <  ? )  =  fi  ( 4 ) 

A  I (3) =At 7) 
ah4»=a{2> 

A I  { E- }  =  A  (  5 ) 

AltS) =A(3> 

A I  (7  )  =  A  t  3  > 

A  I  <  -3  )  =  A  i  6 ) 

AI (9) =A( 9) 

RETURN  -  • 
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..SUBROUTINE  DESUB ( TIME ,X,DX) 

REAL  MACH,KE, IX,IY,IZ,IZX,MASS,K,MX,NY,MZ 
INTEGER  SX,7X 
DIMENSION  X<l>,OXti) 

COMMON  /INT /•  R  J 4) 

COMMON/)PXX/VM2 
COMMON  /TRANS/  A (9) 

COHMON  /OISPL/  U,  V,W*OELXE,!iELYE,DElZ£*XT£tYTE»ZTE»RH 
COMMON  /ANG/fiLFA,BETA,SIGY,SIGZ5DPl,aY2,9P3,OY4 
COMMON  /*A  IR/  MACH , VH  <  VS,  QS,8SQ,QSR,RH0 
COMMON  /FRMO  /  FX,FY,F7,HX,MY,MZ,IX,IZ,IY,IZX,HASS 
COMMON  /OTHERi /IT ERA , OT«  OTMIN, EHAX,SMAX«  TMAX, PRNTI 
CONMON/OTHER2/  THTA, PSI, PHI* KE ,S , 0,SX, IX , K,S, MAXPT 
CALL  EAT08F<A,X(4>,X<9)tX<6J  ,X<7)  I 
UA=V(8> 

VA=X(9) 

HA=X<10) 

CALL  WIND (TIME  »X, OX, UA, VA,WA) 

_CALL  MULT (U,  V,  M  ,  A  ,UA',  VA,  MA) 

Q  «  ’i  i  i  ,  ,  ,  ,  *  *  *  »  *  *  *  »  *  *  *  ♦  *  *  «*'*/.-.'**«'***** 

C  *  COMPUTE  LOSiTARGET  MOTION  USING  EARTH  COORDINATE  SYSTEM  * 

£«*****»***«*»**9»**»**  ?*»*****♦**** 

CALL  TARGET<TXM£,X,D«.' 

OELXE=XTE**X(  11 ) 

OELYE=YTE-XII2 > 

OELZ£=ZTE-X(l3> 

RM=SQRT<OELXE*OELXE*DELY£*DfiLY£«-DELZE*OEL2E) 

SIGZ=ATAN2t3ELYE,0ELXE) 

SIGY=ATAN2 (-OELZE , OELXE) 

C 

£**•***«*«*{-»***•*»***»•**********«» 

C  *  COMPUTE  SEEKER  OYNAHICS  AND  800Y  FIXED  MING  COMMANDS  * 

£#******•*****♦»♦*#***♦****«►»»*****# 
CALL  SEEK£R(TIHE,X,OX> 

£***************»#»»*»»**•**»*««*•** 
C  *  COMPUTE  RHO  ANO  SPEED  OF  SOUND  AS  A  FUNTION  OF  ALTITUDE  * 

C*********»****************»**»***** 
*H=-X (13) 

CALL  AIRTAB(H,RHO,VS) 

VH2=U*U+V*V*W*W 

VM=S9RTCVH2) 

HACH=VM/VS 

C 

V  ♦**********»••*•*•*•  *■***•#■**♦**** 

C  *  COMPUTE  ANGLE  OF  ATTACK  * 

£»***************»**♦*******♦*»»**** 
IFIU.EQ.O.)  U=l.E-?0 
ALFA=ATA-N2!M,U) 

8ETA=ATAN2(V,J) 

(;*»**»******•♦***»»**#**  ************* 
C  ♦  COMPUTE  FORCES  AND  MOMENTS  IN  BODY  FIXEO  SYSTEM  * 

CALL  FOROM(TIME,X,OX) 

C 

C*********************************** 
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#«*=HX-X  f  2>*K< 3>*  {IZ-I':>  -Xt2>  *X<i)*IZX 
AflHY=HY-X(l)*X<  3)*(2X-I.Z>-'(X(1)*X  { 1)-X (3> *X <3>  >*IZX 

- CtFZ=tt2-X-UT»X{?F*-tIY-IX»-X<2T*X<3F*IZX-  . -  ~  ■ 

GEN=IX*IZ-IZX*IZX 

-e  &f>  . . .  . . . . . 

OX(i) =(DKX*IZ+ DHZ*IZX)/DEN 

C  DQ  -  -  -  -  -  --  ---- 

0X(2 ) =CHY/IY 

-  C  OR  -  -—  . 

OX (3 >  =(OHZ*I X*OHX*IZX) /OEN 

— 5  OOHSFRftf-NT- - - — 

KE=tl.-XC4}*Xt4>-Xt5)*X(5>-X(6}*X(6»-xm  *Xt7)  >*K 

€  CEO  *  .  ~ 

0XC4>=-.5MX(5)*X<i>+X<&)*X(2>*X  t  7J  *X  (3>M-t<E*X  (4) 

C  OE1  -••-••  '  '  ~  -  -  -- 

0X(5)=.5*(X(4)*X(1)  +  XC6)*X(3)  -X(7>*X{2>) FKE*XC5> 

€~&€2 —  ' - 

OX(S»  =.5*  (X  {  4)  *X(2>4-X(7)*XC1)-X{5)*X{3)1  +  K£*X  C6) 

C  DEJ  .... 

0X17  )  =.5*  <X(  4>  *X(3)*-X(5)*X{2)-Xl6>*Xtl))  *KE*X  C7I 
c  • 


-  -Ct+Aff&E  -FORCES  FeRrt-BOeY-ftXfX-TO-eiWTH - - -  -  — 


CALL  9FT8EM4,  X<4)  ,*t5f,Xt&»  ,-XT?U  —  -  --- 

CALL  MULT(FXE»FYE,FZe.A,FX,-YtFZ> 

C  ZE  D0UL9LE  DOT  . 

0X<13>=FZE/MASS*G 

— 6-  ZE-  -GOT - - — — - 

0X(13)=X (1C) 

C  ye  DOUBLE  O&T . - . -  - 

OX  (9)  =FYE/M"iSS 

-C  YE-  QOT  -  -  — - - -  -  -  - -  - 

0X<12)=..<9) 

•  -G--XE— OOUBL-E  -BBT - 

0XC8)=FXE/MAS5 

C  X€  flOT  -  --------- . . 

OXlll)=X (3) 

RM=S<lW-F0ELXE*-OELX6*l>EL-VE*9SL-¥€-*-i>£LZ€-«-OEL-ZEJ - -  - 

RETURN 

_  ri]r _ _  - _ _ 

tnir 
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SUfiRGUTINF  RUH-&*(TIME,-y,F,VM)EL,-0£bMI-N,£HAX,-^IG,5HAX,DESO9> 
DIMENSION  VS  <5  0)  ,C<4 ,50) 

01  HE  NS  ION  V(l>  ,F( 1), *(4) 

IFtIC.GT.t  )  30  TO  2  J 

13  H=OFL 

I  FG=  0 
W(1)=H 
W  (2)  =  DEL 
NH=Q 
W  (3)  - NH 
GO  TO  135 

C  REST  ART  IF  INPUT  OFL  HAS  CHANGES 

20  IF40tL.NF.W<  2)  L-  GO  TO  13 

135  TIMEO=TIMF 

0 G  13o  1  =  1, N 

1 J  6  VS(I>=V(I) 

IFdPG.GT.'-)  GO  TO  120 
W  (4)  =1.1 

CALL  0£SUP(TI'1F»V,F) 
no  no  J=  l , n 
1 1 H  C<1,  J)  =r  <J)*H 

12  J  T IN- =  TIN EC+H /2  . 

no  13  j  . J  =  1 ,  N 
133  V  ( J)  =  VSt  J )  +C  (1  ,  J)  /2» 

.4(4)  =2.1 

CALL  OESUO  (T  Ii F. , V  ,F) 

DO  14  1  J  =  1 »  N 
Ifj  CC2,  J)=F(J)*H 
no  15  )  J  =  1 ,  N 
1  j  C  V(J)  =VS<J)0<2,J) /2. 

W  (4)  =3.1 

CALL  DES’JsJ  (THFiViF) 

00  >3  v,  J  =  1 ,  N 
3J3  C(3, J)=" <J)*h 

timi  =  ttm.:o+h 
nc-  3  J  1  J=1 

331  V ( J) =VS( J) +C (3 ,J) 

W { 4 ) =4.1 

CALL  DtSUB  ( T  I'l  E  , V  , F) 

00  5. i  2  J=1,N 

3J 2  C  (4, J)=F (J)*H 
no  3  J  3  J  =  1 ,  N 

533  VIJ)  =VS(J)  MC(  1,J)+2.»C(2,J)  +2  .*C  <3,J)*0  (4,J)  )/6, 

W  (4)  =1.1 

CALL  3C3UP  (T  IMF  »V  ,F) 

00  233  J=1,N 

2 j  i  C!1J)--F(J)’H 

IFC-=1 
IC=1 
RETJRN 
END 
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SUBROUTINE  HUuTIOUTX* OUTY ,3JTZ,3  » INX.INY  *INZ> 
REAL  ENX  » 1NV  »I NZ 
01  HZ  NS  ION  R(9» 

0UTX  =  3(l)*lNXF8*4>*lNY  +  3<7)  *INZ 
OUTl  =  9(2)  *IN:<F'3{5)*INY+‘M8)*IN? 

0UTZ=9(3)-»I  JX4-B(6)  *INY  +  3t'R)  »IN7 

RETURN 

END 


10 


13 


-SUBROUTINE  -RK-t£RtrIM£TVrF»'R«'0£L-*0£LttIN,EMAX.HTfeT3MAX  .DE3U0) 
DIME NSION  V(l)  .H(4> ,C(5.5G) .VS (50 ) »F(i) ,E (SO) 

COMPUTE  STOP  FOR  THIS  INTEGRATION  STEP 
TSTOP=TIMF-OEL 

CHECK  FOR  FIRST  TIME  INTO  ROUTINE  . . 

IFdC  .GT.  0)  GO  TO  20 

€MIN=EM/H(/3eT -  -  -  -  -  - 

FIRST  TIME  IN 

H  =  DEL  ‘  ‘  ‘ 

IFG=0 

NH  =  a 

W  (21  =  DEL 
-GO- TO  105— 


2C 


C 

■C. 

105 


106 

C 

110 


119 

120 


-+30  - 


RESTART  IF  INPUT  DEL  HAS  CHANGED 
IF (DEL  .NE.  RI2TF  GO  TO  13 
H  =  H  (1) 

NH=  M  (31 

INTEGFAT  USING  R-K 

-  -  S-AVE-  -V-  T  A3  IE  *AtrU&S-IN-VS  -ARRA» - 

TIHEO=TIME 
OO  106  1=1. N 

VS(II  =  V(II 
IFdFG.GT.O)  GO  TO  12C 

FIRST  FASS  THRU  R-K  COMPUTATION 

*H4T-=t.i -  •*  '  — - - 

CALL  OESUB  (TIHE.V.FI 

ne-m  j*i.n  -  — 

C(l.  Jl=c(Ji*M/3. 

TIME*TIMEO*H/3 . 

"O  131  J=1 *N 

V  (UF«  VS4U1+G  U.-U - - 

H(41  =2.1 

CALL  OESUB (TIME.V. FI 
00  141  J=l»  N 


140 

150 


310 


C<2tvM=FFU1jl+V-3. - 

00  150  J=1.N 

V  (>!!•- 
M(4)  =3.1 

GALL  OESUB4TIME.V.B  • 

00  301  J=1.N 

C(3,U1=F(J»,»1V3«- - 

TIME  =  TIMFO*H/2  . 

00  3-9 1 — J*hA-  -  — 


301 


302 


313 


V(Jl =VS(J1 +3.*C(1»J1 /8.+9.*:(3»J) /«. 

W (4)  =4.1  .  -  - 

CALL  nESUe(TIME.V.F) 

-00  312  U=1tN  -  - 

C(4,J»=F(J!*H/3. 

•  -TWS-TIMeO+H-  - - -  •  * . 

00  303  J  =  1.N 

V  {  J)  =VS(U)  ♦!  .5*C(l»vl*-4.5*G(  3.  J1  ♦  6.*G(4t-,> 
W(4)  =5.1 

CALt  3ESUP.(TMErVrf>  - - 

00  304  J=i.N 


IFOEL.LE.DELMINl  SO  to  163 
IF(.  5  *H  •  LT  .OE_HIN)  GO  TO  13C 
r  TEST  FOR  HALVING  H 

SMA<  =  ). 

00  153  J  =  1  ♦  N 

:<J)  =(C(ltJ)  -*♦  .3^F(3»J)+4.»:(4fJ)-.5*C(5,J))/5. 
ER»0^=A3S(F( J)  ) 

I F  (ERROR  .  GT.  E'l  AX )  GO  TO  155 
T  f  ( E  R  ROR  .  G  T .  SH  A  X )  SMAX=ERRJR 
153  CONTI  NOE 

IF(5HAX.I.T  .ENIN)  10  =  10+1 
I F  ( S  '1  AX  *  G  h  .  E  MI  N)  IC  =  1 

0  OHEG <  FOR  THIRO  ITERATION  JITH  ERROR  LT  EMIN 

IF{IC.LE.3)  GO  TO  16E 
H  2-Z  .  *H 
IC=1 

I F  (H  2  .  GT  ,  Hf-  l  )  GO  TO  1GG 
H  =  H2 

r  PINISH  TUI?  ITERATION 

GO  TO  15..' 

C  HALVE  »» 

1  5 H  =  t  t*"'  H 

NM  =  NJH+-1 

IC  =  1 

C  START  k'-K  1 4  TE  ORATION  OVER 

T I nE  =  T I M - L 
JO  1  5  7  J=l,  4 
157  V(J>  -  VS(J) 

GC  TO  11  : 

C  iPJATf  V  T  A  3LE 

15  '.  OC  L  3  .i  J=lt  N 

13  V  ( J)  =VG<  J)  t.  ‘J*  (C<  ItJ)  +C(5tJ)  >  *2,  *C(**t  j) 

W  (4)  =1.1 

CALL  OfS'J"  (T  IMF,y  ,F) 

OC  2  1)  J  =  l»  4 
?)'  C  Cl,  J)=r  <J)*H/3. 

IFG=  1 

C  E'lO  OF  REOUt°EO  INTEGRATION 

TC..  W(1)=H 

w(2)  -  DEL 

W(3)  =  NH 

RF  TJ  R  '4 
“NJ 
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c. 

313 


C 

300 


305 


1 

~e 

C 

c 

4 

7 

8 
e 
c 


c 

9 

10 
c 
c 
11 


12 


-strnotm  me--kam-pc<-x — fy-vcert -toeu 
DIMENSION  Y(l)  ,0ERY<1>  «AUX(16,20)  tH(15 
COMMON  IHtF,  ISTEP  - 
CHECK  FOR  FIRST  TIME  .  IN  ROUTINE 
IFtTC.GT.C)  SO  TO  300 
W  (2)  =  H 

IHLF=0 
TIME=X 
GO  TO  305 

RESTART  IF  IN»UF  DEL  HAS  CHANGEfr 
IF  (H  »NE . H( 2)  )  '  GO  TO  313 

GO  TO  211 

00  1  I=l,NOIM  - - - 

AUX( 16*1 ) =0. 

AUX(15,I1=0ERY <I> 

ADXf  1  » I)  =Y  m 


COMPUTATION  0-  OERY  FOR  STARTING  VALUES 

CALL  FCT(X,Y,OERY» 

00  G  1=1, NOIM 
AUX(8,I)=0ERY(I) 

COMPUTATION  Oc  AUXT2.I) 

ISM*1  -  - - -  - 

GO  TO  100 

X=X*H 

B&-*r9  I^ItHOIM- - — - 

AUX(2,I»=Y  (I) 

INCREMENT  H  IS  TESTED  BY  BISECTION 

THLF=T«tF-*-l - - -  - 

X  =  X-H 

"9-ta-  £*trN3  W - 

AUX(4,I) = AUX (2*1) 

-  - -  --  --  - 

N=i 

ISW=2  -  ------ 

GO  TO  100 


13  X=X«-H 

-GALL-  FGT{XtV,4ERY4  - . 

N=2 

BO  1*4  I  =-4, NO  IN  . - 

A'JXl  2  » I)  =Y  (I ) 

14  AUX{9»«-=0€RYJ-T4 -  -  -  -  -  •- 

ISW=3 

GO  TO  180  -  -  -  -  -  . — 

C 

C-  COMPUTATION  OF  TEST  VALUE  OEL-T 

15  OELT  =  0 • 

-OB  tB-I=4rNBW~  - -  - - 
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- OEL-T-OEL  U~MYc4ri$y  - 

OELT  =  .  066666  67  *OELT 
SHAX-  =  OEL  T 

IF  J9ELT-EMSX )  19, 19, 1 7 
i  7  JFHHUF-IOV  11,18,18 
C 

■S  --WS4  TiSCASTOR-Y  --A6CUR-AG-Y--AFTE-R- 14-  -BISECTIONS.  -ERftSR-KSCSACg-.  - 
13  IHLF=11 

-  -  X— X+-^T-  -  ■  -  -  -  ------ 

GO  TO  4 

C  THESE  IS  SATISFACTORY  ACCURACY  AFTER  LESS  THAN  11  BISECTIONS. 

-49 - - : - 

CALL  FCT (X,Y ,OERY ) 

- - 50  24-1^4  ,-NQP< - 

AUX( 3,1) =  Y<I) 

-20 - ftUXt  l- 0,1  CERT- (-T) - 

N=3 

GO  TO  10  0 

C  . -  - 

21  N=1 

X=X*H  ..... 

CALL  ^CT (X,Y ,OERY) 

00  22  I=1,N0IH 
AUX<11,I)=0ERY  <I) 

22  Y ( I) = AUX  <1»I)*N*( .37 5* AUX 18,  I )  7 916667*AUX(9 tl) 
*-.2033333**0X1  U,  I)  ♦  .  04166&67*0£RY  Cl) ) 

23  X=X*H 

CALL  FCT (X»Y ,OtRY) 

-VU2)  =H  -  -  - ...  - -  - 

Wtl)  ~  IHLF 

IC=2  -  ....  . . 

RETURN 

24  .IF<N-4>  2S,2Jd-,204i  - —  - - - 

25  CO  26  1=1, NDIrt 

AUXIN, I)  =  Y  ( U  .  .  ....  . 

26  AUX( N*7, I) -DER Y i I ) 

IFtN-3)  27,29,200 

C 

27  DO  2-a_X=l,U0Irt-  —  --  - - -  _ - - 

OELT  =  AUX  (9,1  )*  AUX  <9,  I) 

CELT  -  DEL  J-*l)£Ll-  --  -  -  -  -  . .  .  - 

28  Y  tl) =«UX (1 ,I>* .  33333  33*HM*JX (8, I)*DELT+AUX{10,I) ) 

DO  TO  23-  . . . . . 

r, 

23  -  OO  3d  - - - 

OELT-  AUX (9,1 )> AUX l 10 . I) 

0"LT=0CLr+OELT*0ELT-  - - - -  -  .  _  ... 

30  Y  (I)  =  AUX  (1 ,1 )  ♦■  .  375*K*  i  AUX  (8, I )  OELT+AUXT  11,1)  ) 

GO  TO  23  -  -  - - -  ... 

i. 

C-  THE  f 3LL OH  M-4-  PART  OC-  $03 ft 00  flNi-  -HAHUC  COMPUTES  -BY  MEANS  OF 
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e - RUNGA-KUTTA-  HE-FHOO-ST  ART  INS  VALJE-S  -F0R~T1)€ - - 

C  PREDICTOR -CORRECT OR  METHOD 

100  00  131  I=1»HBIM 

Z=K*AUX(Nt75 1) 

AUXID,!)-?  -  . 

101  Y(I)=AUX(N,IM-.4*Z 

C  Z  IS  AN- AUXILIARY  STORAGE  L-OO-A-TI ON  -  --  -  -  - 

C 

Z=X+.4*H 

CALL  FCT(Z,Y,0€RY> 

00  132  I  =  1-,NDIM 
Z=H*OERY (I) 

-  -  -  -AOX+6-rI>  «? - - - - 

102  Y  CI»  =  AUX  (N,I)*.2969776*AUXC3,I>«-.i587595*Z 

Z=X+.  455  7372  *P 

CALL  ECT  (Z,Y,3.RY>  -  •  -  -  .  _ 

DO  103  1=1, NOI) 

- 2E»*<v&ERY-H') - - - - 

AUX(7,I)=Z 

103  Y(D=AOX(N,I)+,2i8l0  94*AUXt3,I>-3  i‘v965*AUX  (b,  I)  0.4328*5»Z 
C 

Z=X*H  ..... 

CAL'.  ECT  (Z,Yt3ERYI 

.  OO-lOV-f^lrNOxK— -  - - - 

10  4  Y(I)  =  AUX  'N,I)f  .174763 3*AUA (5, I>-.5514u-.'* AUX .*6,1) 

-$+1.295536»A1JX<7,mr17il84«*«*Qc*YfI>  . —  -  - 

GO  TO  (9,13,15,21)  ,I$H 
200  IST£P=3  - 

231  IFCN-8)  2 04  ,202,204 

2  )2r  -0&-2-0^-NN2t7 - - 

00  203  I  =  1,N0IM 

AUX(  N-l,  T)=kUX(N,II  ..  .. 

203  AUXi  N*-6»  I)  =4UX  (N*-7,I« 

H *7  -  — - - -  -  -  -  .....  .... 

234  N=NU 

DO-245 -I-WWOTM - - 

AUX(M-1,I)=Y  (I  > 

205  AUX(N+f.,I)=OE4Y(I5  -  - . _  .  . . 

X=X+i« 

2t>6  ISTEP= tSYER+l  ...  . _ 

DO  207  1=1  ,ND£  N 

OELT-AOVTN-4^1  H  i-  ;3-i333--*»  I  i 3— AO X  <N*S ,  34--^ 

$AUX(N*4, I) ♦AUX (N+4,I) * 

-  T<I)  *OELT-,3*SDi98*AOXfi6,I> -  -  -  -  - . . . . . 

237  AUX{ 16,1 ) =OEl T 

CALL  FCT  (XrY-,9ERY)  -  -  - . .  - . 

00  2  38  1  =  1  ,MCIM 

OELT  =>-.-l2E*4R  r*  AOX+N-t-c  D— AIK  («->•»  fOER^M ifM-AOX(-N*&d4-f — 

■SAUXCNiStD-AUX  (N*5,I) )) 

A 07  ?  1  6rl  >-=£0  X-L1&,  I>— CELT  - . .  - . 

208  Y ( D =0EL T  +  .O  ?4 330 17* AUX *16,1) 

OFLT  =  0*  •  .  _ _ _  .. 

00  209  1=1, NOTH 

2S9  ouLT  =U€LT+AtiX {•)•?■,  If*-AOSfAOX-{  L6-»1H -  - -  — 
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- 4444E44-E-HA* - > 244-r422  ,422 - — - - 

210  CALL  FCT ( X  ,Y ,3  ERY ) 

IC=1  -  - - -  - -  -  ~~  -  -  -  - 

W<2>  =H 
W  Ml  — M.F 
RETURN 

444 -  ■  -IF<4-HLF-411-  -2l-5r412-r414 - - - - 

212  W(1>=IHLF 

-  --P.E-T4R4  - -  --  -  -  ------ 

215  IF(0FLT-.G2*E*1AX  )  216,213*201 


217  IF1N-7J  201,218,218 


ISTEP=G 

. OC  -2-24-1-=  l  .MOI  M - - - 

AUX( N-l, I ) =AU< (N-2,1) 

AUX{f4-2,  n=AUXUl--4,I4  -  . 

AUX(N-3,I)=AUX  <N-6,I> 

AUX(Ut6,I)=AUX  IN+5.I1 
AUX(N*5,I)=rtUX (N+3,I) 

... _ _ fluy  (  m  j )  -ftL'X  <N*  1 1 1 ) _ _ _ 

0ELT  =  AUX(N40,I)>AUX{N+5,I) 

DE4T-OE1.T  +DEcT  -*-OELT  _  - 

22  i  AUX 116, 11=6-.  95  2963*  (V  (II-AUX  (M-3  »IJ)-3»3&lill  *H*  tOERY  (II  4-DELT 

4*AUX(N*4,U» 

GO  1  0  201 

ZZZ--  IHLF=1HL41*1-  -  -  ~  - . -- .  -  - - - - 

IF  f I HLF- 1 G )  223,223,  210 

223  N  =  -  ..... 

I 3TEP=0 

- CO.  Z2'+  1  =  1  .KOI  M  .  ...  . .  -  -  -  -  -  . — 

Y  ( I)  =  .  00  3962  5*  C  30  . *AL'X  <N-1, 1 1  ♦  13  5  .*AUX  fN-2 , 1)  «-40. *AUX  CR-3 , II  ♦ 

3AUX(M-ri.  ,1)  >  ^.1171875  Y-tAUXLtii6,  IJ  -6,  •AUXlN-t.5,1  J--Ai1X_lH*-4,Il  UtU — _ 

AUX<N-<*,I)  =  .0&39  062  5*U2.*AJXW-1,I}*135.*AUXCN-2,I1  + 

3 1 0  8.  *  aOX  ( N-3  ,1 1  *AUX1  U.-4,  1 14  - . 023  4375*  UUX  l  N*5  18.*  AUX 1 

1  9.*AUY(N  +  4,  II  >*H 
A1!X(  N-3, 1)  =AUX  (N--2.I)  - 

22  4  »UX(N+4, I) =AUX (N*5,I1 

OELT  =X-(H*K) 

CALL  FCT  (DEL  T.-Y..UERY  l-  -  -  . .  . . 

OO  225  1=1, NOTH 

-  -VJX(  N-2, 14=Y  (I) - .  —  .  — - -  - - 

AUXCN  +  5,  I )  =0ER  Y( I ) 

225  -  -Y4IJ  =AUXJLV-vj..I-i - - - 

3£LT  =  0ELT-{H-H  i 

CALL  -ECT  (DDL 1»-Y .QERY4-  .  .  .  -  -  -  - - -  -  - - 

00  226  I=1,K9IM 

OELI=  AIJX  (Ni5-.I1  *A4!X  ( <4*4  ,  I J  - - -  - - _ 

0ELT=  3ELT  *0f LT  *OELY 

.  AWC(  10rl-  -8.  962963M40744*-l-rI->-^4-t4M--3^341.141-li;*<-AL,XlA(»6^I44tO£4D  - 


77 


- SUBROUTINE-  T*Tfc£ -TNrANSHEf  -  T-.-v .  »  TrHX-vNPX-fY  ,4  T-fHY-,  HP  YtZ  ,ZT  ,*Z  i NPZ3- 


C 

r 

c 

c 

c 

-fr 

c 

c- 

r 

d- 

c 

'  </“ 

c 

c 

c 

c 

c 

-0 

c 

c 

c 

c 


c 

o 

c 

c 


53 

52 

51 


45 


42 


3  INDEPENDENT  VARIA8LE5- 
'.Ai  Im31.ES  fOROER) 


TA8LE  LOOK-UP  ROUTINE  PGR  1,^ 

N  =  NUMBER  OK  INDEPENDENT 

N  =1  F 1RST  ORDER  (X)  -  ■  -  ■  ■ 

N  =2  SECOND  ORDER  <X.r > 

-;4  ...  ORDER-  (X,Yr-Z  )  - -  -  - -  - 

ANSWER2  (DEPENDENT  VARIA5LE  CORRESPONDING  TO  INPUTS  X,»,7) 

NT  *  MOLE  OP  DEPENDENT  VARMALE  CORRESPONDING  ■--- 

HT(I.J,:<)  INCREMENT  SUBSCRIPTS  LEFT  TO  RIGHT  WHEN  LGAOING 

X  »- THE  ARGUMENT  OR  INDEPENDENT  VARIABLE  X-  -  -  - - 

XT  =  TABLE  OF  T.NOEP,  X  VALJES  CHUST  BE  IN  INCREASING  ORDER) 

-NX - — •  -NUMBER-  CP  -POINTS'  TH-  -XT- - - - 

NPX  r  NUMBER  OF  POINTS  fO  USF  FOR  X  INTERPOLATION 

¥4  •  »  THE  ARGUMENT  OR  INDEPENDENT  VARIABLE  Y-  - . -  - 

YT  =  1ABLE  CF  INDEP.  Y  VALUES  (MUST  3E  IN  INCREASING  ORDER) 

NV  =  NUMBER  OF  POINTS  IN  YT 

NPY  =  NUMBER  CF  POINTS  TO  USE  FOR  Y  INTERPOLATION 

2--—  a  THE-  ARGUMENT  -09- 1 N&E PENDENT  -VARIABLE — 2- - 

ZT  =  TABLE  OF  INOEP.  Z  VALUES  (MUST  BE  IH  INCREASING  ORDER) 

NZ  =  NUMDER  OF  POINTS  IN  ZT 

NPZ  =  NUMBER  OF  POINTS  TO  USE  FOR  Z  INTERPOLATION 

REMARK  1.  THIS  SUBROUTINE  WILL  ACCEPT  1ST,  2ND.  OR  3RD  ORDER 

--REMARK  -2-.—  IF-  1ST  ORDER-* -USE— XT-  AND-  NT* - 

REMARK  3 «  IF  ZNO  ORDER,  USE  XT , YT  ANO  NT. 

REMARK  4.  F  320  ORDER,  USE  XT  tYT.ZT  AND  NT,  - 

REMARK  5.  ALWAYS  USE  NT (  »  ,  )  FOR  THE  TABLE  OF  DEPENDENT  VALUES. 

DIMENSION  XT  (MX)  ,  YT(N>  )  » ZT ( YZ )  «WT (NX.NY  »NZ)  ,NU5)  ,A(iO) 

GO  Y-G— t5T-r52 ,5-3  -N  •  - -  - - 

CALL  LI^IT  (Z,ZT,NZ,NPZ,KIMZ,HAXZ) 

CALL-  LIMIT  (-Y,  Y-T.  HY,  NPY-,  MINT -.MAC  Y  )  .  —  -  « 

CALL  LIMIT  (X, XT, NX, NPX. MINX, MAX X) 

2=1 

-GtY-TO  N -  -  -  •  — - 

00  41  J=NTfcZ  .MA.tZ 

•  DO  4-2  “I=-NiNY  *N-AXV  •  ....  .... 

CALL  INTEFP  (2  ,X,XT,WTtl,I,.i)  ,NK  ,  NPX  .MINX  ,  MAX  X,H  (  I)  ) 

4NSWER=N{I»  -  -  -  -•  -  -  - 

IF(N.’£Q.  1)  GO  TO 


CALL  INTERP  (2  ,Y,YT,H,Nr,NPY,MlNY,MAXYTMJ>> 
ANSNER=A-WY  --  --  -  -  -  - 

IF(N.CQ. 2)  GQ  TO  4S 


41  CONTINUE- .  — -  -  -  -  -  - 

CALL  INTERp  (?, Z,ZT, A, NZ, NPZ, MINI, KAXZ, ANSWER) 

>.?  ■-  RETURN -  — .  —  -  -  — -  - 

&N0 
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c 

c- 


25 
\  c 

\* 

\  — 


c 

21 


26 

C 

22 


SUBROUTINE -iriN^F-(-Xt  XT  *Ny*Ni>  ♦  HINX-*HAXX1 - 

THIS  SUBROUTINE  MILL  FIND  He  MINIMUM  AND  MAXIMUM  SUBSCRIPTS 
{■OR  RANGE)  T©  BE- CONSIDERED  FORINTERPOtATfON. 

DIMENSION  XTIMX) 

NPX=-HP-  -  -  — -  -  -  ~  - 

IFINPX.GT.NX)  NPX=NX 

-  ©0  -25-  IdyNY - - -  — - 

IFIXTIII-X)  25  *  22,21 

CONTINUE-  -  -  -  -  --  . 

.  GREATER  THAN  N&X  SUBSCRIPT 


MAXX*NX~ 

MINX  =  NX-NPX*1 

.  WITHIN  RAN6E 

MI  NX -  I-NPX/2 —  -  -  -  -  - 

MAXX  =  MI.NX*NPX-1 
IF INAXX. GT .NX)  GO  TO  24 
IFIMINX.GE.l)  GO  TO  26 

,-rr:r»  n  rrn  .mTro-.T .Tmrrr-.m - LESS  THAN  SUBSCRIPT - 

NINX  =  1 

MAXX=NPX 

RETURN 

. .  NO  INTERP  NECESSARY 

MINX=T 

■HA-KM-  - -  -  . . . . . .  -  - 

RETURN 

ENG 


C 

C 

C 

u 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


-  — SUBROUTINE  INTERP  tt:*T,XT*TcYf  ,WX*lfPX?*W**H**XfY»  - 
THIS  SUBROUT  INF  PERFORMS  A  SIMPLE  INTERPOLATION. 


INPUTS 

LMT  =  1  PPOGRAM  WILL  OETERMINE  SUBSCRIPT  RANGE  l MINX  TO  MAXX). 

LMT  =  2  PROS RAH  ASSUNS  THAT  MINX  AND  MAXX  SUBSCRIPTS  ARE  KNOWN. 

X  APGONrNT  OR  INDEPEN3ENT -V  ARIABLC  F6R-  MHICH  ANSV;R  IY)  WILL 

BE  DETERMINED. 

XT  TABLE  OF  INDEPENDENT  VARIABLES. 

YT  TABLE  OF  DEPENDENT  VARIABLES  CORRESPONDING  TO  XT. 

NX  NUMBER  OF  POINTS  IN  XT.YT 

NPX  NUMBER  OF  POINTS  TO  BE  JSED  FOR  INTERPOLATION . 

MINK  MINE-MUM  -K-T  SUBSCRIPT  USfrB  FOR  INTERPOLATION. 

MAXX  MAXIMUM  XT  SUBSCRIPT  USED  FOR  ;  NTERPOLATION. 

OUTPUT 


Y  ANSWER  OR  DEPENDENT  VARIABLE  CORRESPONDING  TO  INPUT  IX). 

DIMENSION  XT (NX) , YTCNX) 

IF CLMT.EQ.2)  GO  TO  13v 


iALt-LlW-Ht,  K-l 


130  Y=YTIHINX> 

IFIMiNX.EO.MAX-X)  GO  TO  106 


Y  =  0. 

70  12)  J*MINX**AXX- 
P=l. 

-  -  ---66-  tl-S  I«*IN»»H»»X - -  - 

IF II . EQ«  J)  GO  -0  110 
i*=P»-<  X-XT  FH-YF-fXT-fJT-xmil 
UO  CONTINUE 

120  -Y cy  HfTFUHM*  . 

100  RETURN 

- EN© - 
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con  o  ooo>  -<■  o  -  oooottifoooooooooooonooorpoooooQo.'^ooo 


SUBRBUTI-ME-QSBSOE - - 

RFAL  MACH , KE  ,< 

INTEGER  SX,TX 
C0MM0N/XXX/VH2 

COMMON  / AIR/  HACH-,VM,  VS,QS.ijSD,QSR,RHO 

COMM  ON/OTHER  2/  THTA,  PSItPHI,  KE  ,S  >  D  ,;X ,  TX  ,K  ,6,  MAXPT 

£  #  ^  4  ^  ■  T .  T  .  _ »  .  «  J*  »  J*  « * — •  _« — *  *  »  *  *  -•  -»  -* 

c  «  commute  dynamic  pressure  terms  * 

c.  ,«•*•«»>***•**»••»  *.*.**»«». **»«*v** 

QS=RH0*VM2*S 

QSC=QS*0 

QSR=RH0*VM*3*3*D 

—  RETURN- -  -  - - - -  --  - 

ENC 


- SUBROUTINE  -RfcO~TNft'vHrttTNl-~ NS) 


PLOT  10 

-  »♦»»-»-»»»-»» wwwwwt wm>w»w»«> »»*»♦»♦♦♦<•«»  2a 

PLOT  33 

-sueRoorrwF -plo-t . .  -  -  - - -  <»o 

'  PLOT  50 

-  -  -  PURPOSE — - - - •  •  ‘  - PLOT — 69- 

PLOT  SEVERAL  CROSS-VARIABLES  VERSUS  A  BASE  VARIABLE  PLOT  TO 

--  - . .  -  - .  PLOT  8C 

U3A3E  PLOT  9S 


CALL  PLOT(A,H,H,NL.NS) 

--OES5RIPT-ION-0P-  PARAMETERS - - - 

A  -  MATRIX  OF  data  TO  BE  PLOTTED)  FIRST  COLUMN  REPRESENTS 
BASE  VARIABLE  AND  SUCCESSIVE  COLUHNS  ARE  THE  CROSS- 
VARIABLES  (MAXIMJM  IS  9). 

N  -  NUH9ER  OF  ROicC  IN  MATRIX— A-  -  - 

M  -  NUN^FR  OF  COLUMNS  IN  HATPIX  »  (EQUAL  TO  THE  TOTAL 

-  -NUMBER  Of  VL-RI  ABLEST  T  -MAXIMUM  -fS-i-Oi - 

ML  -  NUMBER  OF  LINES  IN  THE  PLOT-  IF  A  IS  SPECIFIED,  50 
LINES  ARE  USED, 

N3  -  CODE  FOR  SORTING  THE  BASE  VARIABLE  DATA  IN  ASCENDING 
ORDER 

•J  SQPTTNG  IS  MOT  NECESSARY  (ALREADY  IN  ASCENDING 


-CTH’i:’; -  - 

i  SORTING  is  necessary. 


REtV  .KS 
NONE 

- SUBROUTINES— AUB 

NONE 


OIMNSION  OJT(  101)  ,«PRtll),ANG(9>  ,A(1> 
2,ANC(9)/1L  ,  M*,iH»,iH.,lMX.lHC,  1H0,1HA, lHB.iHC/ 


2  FORHATUh'  ,F11  .4,SX,lGlAi) 

3- FORMAT (1H —  -  -  - 

4  FORMAT  (10H  J‘3LV6/89) 

—  5  -F-ORMA-T-fi-S-'  1 - * - - - *  -  —  —  — - - — 

FGRMAT(.'/3aX , ’HORIZONTAL  AXIS  RESOLUTION  ♦OR  -  •••Fll.A/ 
X  30  X,  *  VERTICAL  AXIS -RESOLUTION  +OR--  »**F11.4l 
FORMS  T{ll»B,liX  ,F12.5,84X,F12.5» 


PLOT  110 
- - PLOT -120 


FLOT 

140 

PI.OT 

150 

PLOT 

160 

PLOT 

IT  b 

PLOT 

180 

*  -  —  PLOT 

199- 

PLOT 

200 

PLOT 

210 

PLOT 

220 

PLOT 

230 

24  0 

-  -  PLOT- 

259- 

PLOT 

260 

PLOT 

2/0 

PLOT 

280 

PIOT 

296 

ploy 

300 

-  -  —PLOT- 

-5*9 

PLOT 

32  3 

PLOT 

330 

♦♦•••PLOT  340 

PLOT 

350 

PLOT 

3/0 

AHSTW- 

PfcOT 

380 

PLOT 

40  0 

•  PuOT 

410 

PLOT 

420 

-  -+4  8T 

459 

NLL=  ML 

IF (NS)  15,  16,  lb 

SORT  BASE  VARIABLE  OATS  IN  ASCEXOINC  ORDER 


80 


o  t>  a  no  o  onaaar>oiv>af> 


- 3r&-&9-44 -IMpt-N- - 

00  14  J=1»M 

. If  tAT44~M«»4-  14~r-14-r -34-  - 

11  L-I'N 

LL=J-N--  -  -  - . 

00  12  K=1,M 

- 4,^4  - 

LL=LL+N 

-  -F=AfL-»- - - -— .  — - 

A(U  =  A(IL.) 

- t2_A4Wt-^=f— - - 


14  CONTINUE 


CONTINUE 


FIN3  SCALE  FOR  BASE  VARIABLE 
-  XSCAL~<A4N1-Al  VDV4fL0A-T  (NL«--3Uf- 
FINO  SCALE  FOR  CROSS-VARIABLES 


H1=N*1 

YNlN  =  AtMt>-  - - - - 

YMAX  =  YMIN 
H2=H*N 

00  4  3  J=H1,N2 
IFIA(J)-YNIH)  2-B, 26,26  - 
26  IF (A { Jl-YMAX)  40.40,30 

24  YHIV=A(JI -  -  — - . - 

GO  TO  40 

30  YMAX=A(J»  -  ... 

46  CONTINUE 

YSCAL=(YNAX-YNINJ  /100.3 

F I  NO  BASE -VARIABLE  RRIJLT  POSITION 

XB=A  41 —  -  - 

t  =  l 

NY=f1-i  —  - 

1=1 

45  -F=I-3r - -  -  - - - - 

xpr=xb+f*xscal 

IF tiULl—XPRl  S4.E0.T3  -  -  -  - 

—  FIND- -CROSS-VARIABLES . - 

- *4-  00-45  INa.UlOl-— - 
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55  -OOT<  FM=flfc*MK 

IF (I. EQ. 1.0%. I .EQ.NLL)  GOTO  U 
OUT(n=OUT(10i  *=ANGt  3  > 

GO  ro  102 

131  00  103  1  =  1,101 ,10 

1 J  0  0U1 ( I)=ANGf3) 

19  2  CONTIWUr 

00  6-3  J=l, HY 
LL=L  +  .J*N 

JP=t  (  AtLL) -YMI  N) /YSCAL3+1.G 
OUT<  JP)  =  ANG<  J) 

GO  CONTINUE 

IF  <1  .  rttr.  NLL>  -  GO  TO  2uc 
IF(Y  MIN.  GF  .0  .)  GO  TO  20U 
IF  (Y  MAX.  Lf  *3  •>  GO  TO  203 
IZfR3=-YHlK/YSC.AL  +  l. 

OUT( I  ZERO) =ANG (5) 

230  CONTINUE 
C 

C  PRINT  LINE  ANO  CLEAR ,  OR  SKIP 

C 

WRITE  (6,2)  X?R,  (OUT ( I Z )  ,IZ  =  1,101) 
L  =  L*  1 
GO  TO  80 
7 1*  WRIT  £  (6,  3 ) 

3>:  1  =  1*  1 

TF(I-NLL)  «*5,8A,86 
P.u  x  PR=  A  ( N) 

GO  TO  59 
C 

C  PRINT  CROSS-VARIABLES  NUMBER'S 

C 

86  CONTINUE 

YPR( 1)=YMIN 
00  90  KN  =  1  ,9 

93  Y PR( <N*1 ) =YPR{ KN) *YSCAL*lU.3 
•YPR(  11)  =  YN  AX 

WRIT  £  (6, 8 )  YH? (1)  ,YPR  (11) 

WRIT E  (6,7)  YS  CAL ,XSCAL 

RETURN 

END 


SDftROUTf-NE  trAGNINPWiOttf^OTVX'T^  i-fNOEX  ifl  MCONt 
REAL  INPUT 
DIMENSION  X(  1)  ,DX  (1) 

DX(INOEX)  =  (-<  (INDEX)  ♦ INPUT)  /  TIM CON 
0  UT3  U  T=X  ( I  NO  £X  ) 

P£TU  ?N 
END 


SUBROUTINE  SEC  ORiH INPUT » OUT3UT  .<  » OX» INDEX  *£T  A  ,WN) 
DIMENSION  X<  1)  ,DXfl) 

REAL  INPUT 
W2=WN*WN 

OX(INUEX)  =X(  INOEX  +  1) 

DX(lNDPX  +  l<  =  -2  .*ETA»WN*X  (J.NOEX+1)  -W2*X  (INDEX)  +W2*  INPUT 
OUT3 UT=X  (INDEX  ) 

RETURN 

END 


SU0R3UT I NE  L IM ( INPUT , OUTPUT,  XMIN , XMAX) 
REAL  INPUT 

IF(TN°UT.GT.  X'lAX)  OUTPUT=XMAX 
IF (I NPUT •  LT.  X1IN)  OUTPUT=XMIN 
PfTURN 
END 


SUBROUTINE  L IN  S  TA ( IN°UT » OUT3UT  tXH IN ♦ XMAX ) 
PFAL  INPUT 

I F ( I N  PUT . G  T . Xi A  X )  OUTPUT  =  IN3UT=<MAX 
IF (INPUT. LT. XT  IN)  OUTPUT  =  IN3UT=X  MIN 
RETURN 
r  NO 


SUBROUTINE  OEA  DSP  t  INPUT  ,  OUT3UT  OwER, UPPER) 
PEAL  INPUT 
REAL  l.OWEP 

Ip (I NPUT .GT. LO HEP. AND , INPUT. L T .UPPER)  GO  TO  i 
TF (I NPUT . Lf. L3HER )  GO  TD  2 
OUTPUTS  NPUT -JPPER 
RETURN 

?  OUT3  UT=I  NPUT -LOWER 

RETURN 

1  OUTPU  T=0 • 

RETURN 

END 


1 


REAL  INPUT 
REAL  MAXLROtMAXNLO 
IFCABSCINPUT).  GT.FOV)  GO  T3  1 
OUTPUT-INPUT*3  LOPE 

IF  (A9S (OUTPUT) .GT.MAXLRO)  OU TPU T=SIGN ( HAXNLO , OUTPUT) 

RETURN  -  •  -  .  - 

OUTPUT=0  . 

RETURN 

END 


SUBROUTINE  3 YR02CT&Y  ,  TG2 , X,  DX  ,  INOEX,  IX  ♦ITP,ITY,WST 
RFAL  LY,LZ 
REAL  IX,ITP,ITY 
DIMENSION  X  ( 1)  ,  )X  (1) 

TMTAG  =  X(INr~X) 

PSIS  =  X ( I NPLX  +1 ) 

CHG=COS ( PSIS ) 

SPG=  S INC  PSIS ) 

CTG=CDS(THTAG> 

STG=  S INC  TH TA G) 

P^X{ i) 

Q=X( 2) 

RP=X-<3) 

LY=TGX*CpG+TGZ  *STG*SPG 
LZ=TGZUTC 

P3=f P*CTG-PR*S TG) /CPG 
QG=-LZ/<PGMIX-ITP)4-W3*IX> 

RG=LY/{PG*(I X- ITYJ  +WS*IX) 

OX  (INDEX)  =OG/3PS+PG*Sf»G-Q* 

DX(INDEX+1)  =  RG-RR*CTG-P*ST3 

RETURN 

END 


SUBROUTINE  LOL  AG t INPUT ,OUTPJ T , Xt OX , INDEX , TMCON1 ,TWC0N2) 
REAL  INPUT 
DIMENSION  Y<  1)  ,QX  (1) 

OUT?UT=(X (INDEX) +THSON1* INP J T ) /TMCON2 
DX  Cl  NDEX  )  =IN  PUT-OUTPUT 
Rfc  TJRN 
END 
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